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ABSTRACT 



The objective of this thesis is to investigate and 
evaluate the effectiveness of modern estimation methods 
with different array geometries as they apply to the 
problem of bearing estimation. The algorithms were 

selected from those that apply to multidimensional case, 
including MUSIC, PHD, minimum norm, and Capon's beam- 
former. These four techniques were chosen based on their 
high resolution capability, and their ability to deal with 
three-dimensional non-uniform arrays and to estimate both 
azimuth and elevation angle of arrival (AOA) . Computer 
simulations were run for linear arrays, circular arrays, 
and combinations of the two. The test conditions included 
(1) two closely spaced emitters, and (2) various levels of 
additive white Gaussian noise. 
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INTRODUCTION 



I . 

A. OVERVIEW 

The angle of arrival (AOA) is a valuable parameter to 
be used in deinterleaving signals because a target can not 
change its position rapidly. Even an airbone radar can not 
significantly change its position in the few milliseconds of 
the waveform pulse repetition interval (PRI) time. As a 
result, the AOA measured by an intercept receiver on the 
radar is relatively stable value. It requires that a large 
number of antenna elements and receivers must be matched, 
either in amplitude or in phase in order to measure the AOA. 

Direction finding (DF) is the area of electronic- 
warfare support whose objective is to find the AOA [Ref. 1] . 
DF systems passively examine the spectrum use by hostile 
emitters and process the signals to obtain enemy AOA 
information . 

The two common approaches to measure AOA are based on 
amplitude and phase comparison [Ref. 1] . The phase 
comparison system usually can generate AOA resolution of ±1 
degree, which satisfies the modern EW requirement. The 
number of antennas that can be used in an AOA system is very 
low in comparison with electronically scanning antennas in a 
radar system. In an airbone system, the maximum number of 
elements might be ten. In a ship-based system, the number 
might be larger. 
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For the phase comparison method at least two separate 
antenna elements with a predetermined space between them are 
required. Since the distances from the emitter to the two 
elements are not the same (except in the broadside case) , 
the incident wave arrives at the two elements after 
traveling uneven paths, and thus it arrives with a different 
phase . The phase delay is proportional to the antenna 
spacing, wavelength and angle of incidence. Since the first 
two factors are constant, the direction information can be 
obtained . 

DF systems often receive signals from several sources 
in the same frequency range at the same time and from many 
different angles of arrival. To resolve multicomponent 
wavefields, modern techniques of spectral estimation are 
used. Most recent digital AOA studies have concentrated on 
high-resolution approaches, such as the MUSIC, ESPRIT and 
minimum norm methods [Refs. 2-4] . In this thesis the 
performance of several modern estimation methods is 
evaluated for different array geometries. This research 
deals with both the one- and two-dimensional problems: 
azimuth angle of arrival for linear arrays, and azimuth and 
elevation angles for a planar array. In addition, an 
underlying assumption is that the receiver is far enough 
from the transmitter and the antenna interactions with the 
platform are controlled so that the arriving wavefront is 
essentially planar. 



2 



The main function of a direction finder is to determine 
the direction of arrival of an incident electromagnetic wave 
relative to the coordinates of the direction finding site 
[Ref. 1] . Figure 1 shows a representative direction finding 
spatial coordinate system with the DF antenna located at the 
origin. The angle of arrival is specified by the azimuth 
angle measured from the x-axis, and the elevation angle 
measured from the x-y plane. Under ideal circumstances the 
electromagnetic field incident on the DF antenna exhibits a 
locally plane wave structure with linear polarization. The 
direction of propagation is indicated by Poynting's 

vector, I 

P = -^Re{ EjxHj*} (1) 



where * denotes complex conjugation, Ej is the incident 

electric field intensity, Hj is the incident magnetic field 
intensity. 



Hj = 



k x Ej 

~zT' 



( 2 ) 



where k is unit vector in the direction of propagation and 
Z 0 is impedance of free space. 

In practice, the incident electromagnetic field is 
usually nonplanar with phase- front distortion created by 
multipath and scattering, and depolarization produced by 
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nonuniform ionospheric propagation effects if such a path 
exists. A full-capacity generic DF should measure angle of 




Figure 1. Direction- finding spatial coordinate system with 

the array located at the origin 

arrival in three-dimensional space. However, many 
operational situations require only measurement of the 
direction component in the azimuth plane. 

A plane wave, incident at an angle other than the 
normal to the baseline between the two elements, arrives at 
one element earlier than the other. Figure 2 illustrates the 
basic phase delay concept . An incident plane wave arrives at 
an angle 9 at antenna 1 inducing a voltage which can be 
expressed in complex notation as Fjexp(_/co/) . After 
propagrating the distance t/sin9 , the incident plane wave 



4 



induces a voltage V 2 exp(/co/ - ct>) in antenna 2, where the phase 
delay given by 

<1> = (2c/7t/l)sin0 . (3 ) 

Therefore the bearing angle is encoded as a function of 
phase delay. For the phase difference technique, the bearing 
angle is computed by using Eq. (3), where phase difference 




Figure 2. Phase-difference DF Parameters 

is measured and the baseline distance d and wavelength A 
are known . 

In the case of linear or planar arrays, a phase 
difference can be measured or computed for each combination 
of elements in the array. Then an array correlation matrix 
can be defined and used in one of several spectral 
estimation methods to estimate the AOA. 
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B. THESIS OUTLINE 

Chapter II summarizes several common antenna array- 
geometries for one-dimensional DF, including uniform and 
non-uniform linear arrays. The combination of circular and 
linear arrays for 2-D DF is also investigated. Chapter III 
covers the mathematical development of the important 
estimation methods, including multiple signal classification 
(MUSIC) , the estimation of signal parameter via rotational 
invariance technique (ESPRIT) , minimum norm, conventional 
and Capon's beamforming, and Pisarenko harmonic 
decomposition (PHD) . Also, DF for multiple emitter AOAs is 
introduced by using the array covariance matrix. Chapter IV 
describes the development of a computer code to simulate 
direction finding antenna performance for signals with white 
noise present . Matlab is used to implement the estimation 
methods and perform the array matrix calculations. A 
graphical user interface (GUI) was developed to allow a 
variety of parameters to be changed with a minimum of effort 
on the part of the user. Chapter V shows the simulation 
results, and validates the code by comparing the data with 
previously published results. Finally, Chapter VI presents a 
summary, conclusion and recommendations for future work. 



6 



II. 



ANTENNA ARRAYS 



Antennas can be viewed as spatial filters; they enhance 
signals in desired directions while simultaneously 
suppressing signals in other directions. In a phased array, 
the selectivity is achieved by phase shifting and then 
superimposing the outputs of all elements [Ref. 1] . For a 
passive array, the improvement in the signal-to-noise ratio 
(SNR) is potentially equal to the number of elements. In an 
array of identical elements, there are two related factors 
that must be considered. The first is the geometrical 
arrangement of the array (linear, circular, etc.) . Second is 
the relative displacement between the elements, once the 
element arrangement is set . 

A. UNIFORM LINEAR ARRAY 

In most EW applications, it is desirable to use as few 
elements as possible. For example, DF arrays on a ship are 
linear arrays placed in the horizontal plane to measure the 
azimuth angle. One linear array can theoretically cover up 
to 180 degrees of azimuth angle, although it is often 
limited to 120 degrees to avoid operating in the end-fire 
mode [Ref. 1] . If the elevation angle is also of interest, a 
linear array in the vertical direction has to be added. As 
far as the AOA measurement is concerned, the two linear 
arrays are usually treated separately. 

Figure 3 shows a linear array with N elements indexed 
from q =0 , 1 , . . . , N-I along the x direction. Let us assume 
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that there is only one plane wave of frequency / and that a 
wavefront arrives as shown. If the input signal is a 
sinusoidal wave, then the q th element will have the output 

g(q,t) = Acos[2-nf{t-T )] 

(4) 

where A is the amplitude of the incident signal, and z q is 
the time delay at the qth element with respect to the 
reference antenna element (q=0). Since the distance between 
the antenna elements is d , the delay time can be written as 

- ad sin 9 

r, = (5) 

c 

where 9 is the incident angle of the input signal, and c is 
the speed of light. 




Figure 3 . A linear array with a plane wave incident 
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The minus sign in Eq. (5) is due to the fact that the 
equiphase plane arrives at the q th element before it 
reaches the reference element. Substituting Eq . (5) into Eq. 

(4) gives 

qd sin0 

g(q,t) = Acos[2fn(t + )] . ( 6 ) 

c 

Let k be a unit vector pointing in the direction of the 
wave propagation and x a unit vector along the direction of 
array axis. Eq. (4) can be written as 

g(q,t) = Acos(2ftn- 2nfqk ■ x—) (7) 

c 

where k ■ x — cos(^- + 0) = -sin0 . 

If there are M signals, the output from the q th 
antenna element is a superposition of all the received 
signals 

g(q.t) = X 4 , */„< - — ' — ) ( a ) 

rn=l C 

where f m and k m are the frequency and direction of the m th 
signal . Often exponential form is used to express this 
result as 

M ^Tlf cjdk * X 

g(q,0 = I A m exp(j(2nf m t - - Jm - } - m )) . ( 9 ) 

m—l C 

In order to obtain digital data, the output is sampled at 
discrete times and therefore the time t will be replaced by 
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integers n =0 , 1 , . . . , b - 1 , where b is the number of time 
sample points. If we assume the frequencies of all signals 

are the same, the factor e j2 ^ m ' can be suppressed, and a 
standard phasor quantity results. 

For an AOA measurement, the elements cannot be spaced 
too far apart or an ambiguity will result [Ref. 1] . The 
shortest distance between two antennas must be less than 
half of the wavelength of the highest frequency in order to 
fulfill the Nyquist theorem [Ref. 1] . Therefore, one must 
sample the incident wavefront twice per cycle. 

B. NON-UNIFORM LINEAR ARRAY (MINIMUM -REDUNDANCY) 

By making use of a theorem by Caratheodory [Ref. 5] it 
is shown that, 'for a given number of elements N , there 
exits a distribution of element positions which results in 
superior spatial spectrum estimators than the uniform linear 
array. In addition to less hardware, certain nonuniform 
arrays are capable of large dynamic range (spectral peak to 
background level), lower sidelobes, and the relatively small 
estimate bias values. A special type of nonuniform array is 
the so-called minimum redundancy array (MRA) , in which 
integer spacings of the base spacing occur only once. A MRA 
is capable of resolving N(N - 1)/2 multiple (simultaneous) 

signals, whereas a conventional array can only resolve N-l . 
Some minimum- redundancy array configurations are shown in 
Table 1. A dot in column 3 represents an element and the 
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integers the number of base spacings between them. (The base 
spacing is typically X/2 . ) 

A disadvantage of minimum-redundancy array is that its 
resolution is not easily increased except by adding more 
antennas and rearranging the array to the optimum 
configuration for the new number of antennas. If there are 
M signals, the output from the q th antenna element is 

m 2nf m d„k„ ■ x 

g(q, t)= I A m exp0(2nfj q - )) 

m=l C 

where d q is the distance between first antenna and the q th 
antenna . 



Table 1 . Some minimum redundancy array configurations . 



Number of ..elements . N 


Array length 


Configuration 


3 


3 


.1.2. 


4 


6 


.1.3.2. 


5 


9 


.1.3. 3. 2. 


6 


13 


.1.5. 3. 2. 2. 


7 


17 


.1.3. 6.2. 3. 2. 


8 


23 


.1.3. 6. 6.2. 3.2. 


9 


29 


.1.3. 6. 6. 6. 2. 3. 2. 


10 


36 


.1.2. 3. 7. 7. 7. 4. 4.1 



C. CIRCULAR ARRAY 

A planar array can be used to measure both the azimuth 
and elevation angles. These types of arrays occupy more 
space than the linear array but they might be necessary for 
some special aircraft and shipboard application. Planar 
arrays include circular and rectangular (of which crossed 
linear arrays can be considered a subset) . A circular 
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antenna array is shown in Figure 4 . The radius of the circle 
is li 0 and there are N elements from q =0 to N-l. The 
array is in the x-y plane and the first antenna element is 
the reference element ( q =0) and it is on the x-axis. The 
location of the elements are at integer multiples of the 
angle 4> 0 = 2n / N . 

Assume that a plane wave with amplitude A, and 
frequency fj is incident on this array with azimuth angle <{>/ 

and 9/ (the angle between the incident ray and the normal of 
the array) . The output induced at q th element can be 
written as 

g(q,t) = A, exp(j(2nfj(t - —))) (11) 

c 




Figure 4 . Circular antenna array 

where the space delay i qI referenced to the center of the 
array is 
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x q , =-R 0 sinG/ cos(J) ?; . 



( 12 ) 



In this case, because the antenna locations are equally- 
spaced in azimuth, the angle § ql can be written as 

=# 0 • (13) 

Therefore 



g(q,t) = A, exp U2nf,(t + - 



R a sin 0/Cos(^- -<(>,) 



-)) 



= A,exv(j2n(f,t 



+ 



/? 0 sin9 ; cos( 

77 



2nq 

~N~ 




(14) 



The incident direction of the input signal can be 
represented by a unit vector k, in Cartesian coordinates as 

kj = -(sin#/ cos^/i + sinG/ sinfjj/y + cosG/z) (15) 

where x, y and f are unit vectors in the Cartesian system. 
The minus sign in front means the vector is pointing to the 
origin. The position vector of the q th antenna element is 

r q where 

f q = cos(q$ 0 )x + sm(q$ 0 )y (16) 

so that 



k, -r q =-sinG/ cos(g<(> 0 —<()/>- (17) 

Extending this to cover multiple signals ( m = l,2 M) gives 
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( 18 ) 



g(q,0= X A m exp(j2n(f m ( - R ° km ^ )) . 

»,=/ k m 

The linear array can only provide azimuth angle 
information because the corresponding propagation vector is 
two-dimensional. The circular array can provide both azimuth 
and elevation information since the corresponding 
propagation vector is three-dimensional. 

D . SPIRAL ARRAY 

We examine a variation of the circular array, Eq. (18), 
where the radial distances of the elements increase in steps 

of R 0 / { N - 1 ) . Therefore R q = R 0 q / N is the distance between 
q th antenna element and the center point. The element 
locations lie on an Archemedian spiral. In a circular array, 
all antenna elements have the same R q . The output of q th 
antenna element is as follows: 

m R„K„ ■ r 

g(q,0= ZA m expO'2n(f m t — q — q -)). (19) 

m=l k m 

E. COMBINATION OF LINEAR AND CIRCULAR ARRAYS 

A three dimensional array can be constructed by 
combining planar and linear arrays. Figure 5 shows a 

circular array in the x-y plane and a linear array along 
the z axis. We examine the special case where the elements 

of the circular array have <j> 0 = 90 ° , which results in two 

elements on the x and y axes equally spaced d from the 
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center. Let there be 3 the elements along the z axis, also 
spaced d . Therefore, the array has seven elements, all 
equally spaced. 




Figure 5. There dimensional seven element array 

The positions of the 3 elements of linear array are 
given by the position vector 

r h = (h - l)dz 

where h = 1,2,3, so that 

k m -r h = -cos(Q m (h-l)d) . ( 20 ) 

If the array receives M signals, the output of /?th 
vertical antenna element of the linear array is 

Sl.nearOl. 0 = Z An CXpOV^fj * " ^ ~ » ( 2 1 > 

m=l K m 
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Another four antenna elements comprise the planar array 



which has an output of the form 



S circular (Q< 0 



S A m exp(j2n(f m t - 

m=l 






)) 



( 22 ) 



where q=l,...,4, M is number of signals, and k m -r q is the 

same as Eq. (18) . This array combines linear and circular 
arrays to provide 2-D direction finding. From Eq. (21), it 
is evident that the vertical array only provides elevation 
angle information. 

We have described the various array geometries . The 
elements comprise the data acquisition system which samples 
the incident signal wavefronts. The next step is to carry 
out the estimation process. In the next chapter, we 
introduce spectral -based algorithmic solutions to estimate 
the AOA. 
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Ill . SPATIAL SPECTRUM ESTIMATION 

Many different high-resolution approaches can be used 
to estimate AOAs from digitized input data. In this chapter, 
five high resolution methods will be discussed: 

(1) Multiple Signal Classification (MUSIC) [Ref. 2], 

(2) Estimation of Signal Parameter via Rotational 
invariance (ESPRIT) [Ref. 3] . 

(3) Capon's beamformer (Maximum Likelihood Method) 

[Ref . 6] . 

(4) Pisarenko Harmonic Decomposition (PHD) [Ref. 7] . 

(5) Minimum norm [Ref.l]. 

All of these are cabable of (1) high resolution, (2) 
estimating both azimuth and elevation AOA, and (3) handling 
three-dimensional non-unif ormly spaced array signals. 

A. MUSIC 

MUSIC is a technique used to determine the parameters 
of multiple wavefronts arriving at an antenna array from 
measurements made on the signals received at the individual 
elements [Ref. 2]. 

The waveforms received at the N array elements are 
linear combinations of the incident wavefronts from M 
narrowband signals and noise. Thus, the multiple signal 
classif ication approach begins with the following model for 
characterizing the receive vector G as 
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G= AF+V 



(23) 



where V a is noise vector. Expanding Eq. (23) 
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The size of each matrix is as follows: 

G : vector of element outputs, TV by 1 

A: matrix of mode vectors, N by M 

F: vector of signals, M by 1 
V : noise vector, N by 1 

The incident signals are represented in amplitude and 
phase at some arbitrary reference point, for instance, the 
origin of the coordinate system, by the complex vector F . 
The elements of G and A are also complex in general. The 
a„(0 m ) are functions of the signal arrival angles and the 

array element locations. That is, a n (Q m ) depends on the nth 

array element, its position relative to the origin of the 
coordinate system, and its response to a signal incident 
from the direction of the m th source. The m th column of A 
can be considered as a 'mode' vector of responses to the 
direction of arrival 0 m of the mth signal. This N by 1 

mode vector will be denoted by a(Q ) . 
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The M by M covariance matrix of the G vector is 

S=E[GG*] = GG' =AFF*A' + VV' (25) 

where E is expectation operator and the asterisk denotes 
complex conjugation. Eq. (25) can be rewritten as 

, A < 

S = APA + XS Qr where P = FF' is a diagonal M by M matrix and 
S 0 is the noise covariance matrix. 



When the number of incident wavefronts M is 
the number of array elements N , then APA' is si 



has a rank less than N ; therefore 



APA 



= S-ASJ = 0 



less than 
ngular. It 
This 



equation is only satisfied with X (not to be confused with 
wavelength) equal to one of the eigenvalues of S in the 
matrix of S 0 . Therefore X can only be the minimum 
eigenvalue X min . Observe that any vector orthogonal to A is 

an eigenvector of S with the eigenvalue a 2 , the noise 
variance. Hence, we can write 

S = UAU h = U S A S U? +U n A n U“ (26) 

where A = diag {Xj,X 2 ,X M } is a diagonal matrix of real 

eigenvalues. The superscript H represents the Hermitian, 
U s is the signal eigenvectors and U N is the noise 

eigenvectors. The noise eigenvectors can be used to form an 
estimator for the spatial spectrum 
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( 27 ) 



Although MUSIC 
of the array 
possible angles 



Pmi(Q) a(Q) fI U N U kCi(Q) 
is highly robust it requires characterization 
response and a spatial search through all 
of arrival , 0 . 



B. ESPRIT 

ESPRIT [Ref. 3] eliminates the undesirable features of 
MUSIC, but only under certain constraints. Assume that there 
are N points of data (N antennas), which can be divided 
into two groups, each consisting of N-l points. The first 
group has an output signal x(t) = As(t) + n,(t) and the second 

group y{t) = ABs(t) + n 2 (t) , where A is a steering vector (not 
related to A in Eq. (24)), s(t ) denotes the baseband 
waveforms, and B is a diagonal matrix 



V; 0 0 

0 v 2 0 0 



B = 



[0 0 0 0 v 



M 



(28) 



where v ; =exp(— j — <7sin0 ( ) and d is antenna spacing. These data 

A 

will be used to make two matrices by the covariance 
approach . 
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Let R xx be the auto-covariance matrix of x(t) and R xy 
the cross-covariance matrix of x(t) and y(t) . If there are M 
signals then 



R xx = E[x{t)x' (t)] = ARA + cr 2 1 (29) 

R xy =£[x(/)/(0] = ARB' A* +cr 2 J 1 (30) 

where £=[$(0*? (0] is the signal covariance matrix under the 
assumption of uncorrelated noise, / is the identity matrix 

and J i is an N by N matrix with ones along the first lower 
diagonal off the major diagonal and zeros elsewhere. For a 
pairwise matched array with doublets the cross covariance 
matrix becomes 

R xy = ARB' A’ . (31) 

In the absence of coherent signals, eliminating a 2 , we 
obtain 

C xx =R xx -<x-/ = ARA 
Cxy=Rxy- cr 2 j, = ARB' A' 

which gives 

Cxx" 7 C X y= AR ( I -yB')A' . (32) 

The singular values of Eq. (32). are given by the roots of 
\l -yB' 1 = 0. The desired singular values are 
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Yk = v * / 



k=l ,2, . . ,M. 



To find the angles we note that for a uniform linear 
array, the angles of the eigenvalues are equal to 2c/7isin0 so 
that 

0* = s in -/ (q-y-) (33) 

2dn 

where 9 k is obtained from y* = exp(jS k ) . Thus, the direction 

of arrival can be obtained without a search technique. In 
this respect computation and storage costs are reduced 
considerably. 

C . BEAMFORMING 

1 . Conventional Beamf ormer 

The first attempt to automatically localize signal 
sources using antenna arrays was through beamforming 
techniques. The idea is to steer the array beam in one 

direction at a time and measure the output power. The 

steering locations which result in maximum power yield the 
DOA estimates. The array response is steered by forming a 
linear combination of the antenna outputs 

y{t)=iw] gl (t) = W H G{t) (34) 

;=/ 

where g,(/) is the output of the / th element, and w,(/) the 

corresponding complex weight factor. The output power for b 
time samples is 
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( 35 ) 



/ b-l. 2 

pm = TZ\y( 0 

b t=o 



= - 1 '}V H G{t)G H {t)W = W"RW 

b t=o 



where R is the signal covariance matrix defined previously. 

The conventional beamformer is a natural extension of 
classical Fourier-based spectral analysis to antenna array 
processing. Suppose we wish to maximize the output power 
from a certain direction 0, given a signal s(t) emanating 
from direction 0 . A measurement at the element output is 
corrupted by additive noise «(/) and written as 



ct.(/) = a i (Q)s j (t) + rtj(t) . The problem of maximizing the output 

power is formulated as, 

max[E{jV H G(tJG H (t)w} ]= mzx[W H E{G{t)G H (t)\v ] 

= max^E\s(t)\ 2 \w H a(Q)\ 2 + <r 2 |wf j ( 3 6 ) 

where the assumption of spatially white noise is used. The 
resulting solution is 



Wbf = 



a(B) 



yja H (Q)a(Q) 



(37) 



The above weight vector can be interpreted as a spatial 
filter, which has been matched to the impinging signal. 

Using the weighting vector, the classical spatial 
spectrum estimator is 

P BF (Q) = a H (Q)Ra(Q) (38) 

where R is the signal covariance matrix. 
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For a uniform linear array of isotropic sensors, the 



steering vector is 



where O = -kd sin 0 
transpose. The 

linear array is 



a(9) = 



1 e* 









(39) 



, N is number of elements and T denotes 
standard beamwidth [Ref. 6] for a uniform 
2n 

— , so sources whose angles are closer than 

N 



this beamwidth will not be resolved by the conventional 
beamformer, regardless of the available data quality. 

2. Capon's Beamformer 

In an attempt to alleviate the limitations of the 
conventional beamformer, such as its resolving power of two 
sources spaced closer than a beamwidth, an improved method 
was proposed by Capon [Ref. 7]. The optimization problem was 

posed as min{P(W )} , subject to W H a(Q) = 1. Hence, Capon's 
beamformer (also known as the minimum variance 
distorsionless response filter) attempts to minimize the 
power contributed by noise and any signals coming from 
directions other than 0, while maintaining a fixed gain in 
the look direction 0 . The optimal W can be found using 
methods such as the technique of Lagrange multipliers [Ref. 
6] , resulting in 



W C AP = 



R~'a(S )) 



^7 // (0)^” / a(0) 



-l 



(40) 
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Using the above weight vector, we can find the following 
spatial spectrum 



^CAp(^) ~ 



1 



(41) 



The power minimization can also be interpreted as 
sacrificing some noise suppression capability for more 
focused rejection in the directions where there are other 
sources present. The spectral leakage from closely spaced 
sources is therefore reduced, though the resolution 
capability of the Capon beamformer is still dependent upon 
the array aperture and on the signal-to-noise ratio though 
R . 



D . PHD 

The PHD is based on the use of a noise subspace eigen- 
vector to estimate the incident angle and assumes that the 
process consists of M complex signals in additive complex 
white noise. It derives the angle of the signals, their 
powers and the white noise variance from the known 
autocorrelation sequence. This algorithm accounts for the 
conjugate symmetry of the eigenvector. This method is based 
on an estimation of angles by using the orthogonality 
property between the signal vectors and vectors in the noise 
subspace. The minimum eigenvector is orthogonal to the 
signal vectors, thus we have 

a, H (9)w ;V = 0 , i = 1,2. M . (42) 
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Once the minimum eigenvector u N is found the from U N in Eq. 

(26) , the angles of the signals may be determined by 
factoring the eigenfilter polynomial [Ref. 7] . The roots are 
guaranteed to be on the unit circle. The a(0) is a steering 

vector that will be zero whenever 0 = 0| , one of the angles 
of the signals. This leads to the angle estimator function 



PHD 



(0) = 



a H (Q)u n 



(43) 



E. MINIMUM NORM METHOD 



This method and the MUSIC method are somewhat alike. 
The basic idea is to find a vector D that is a linear 
combination of eigenvectors in the noise subspace. This 
vector D can be written as 



D-{ &0 6, . . . 8„ ] T = (l 5, . . 5 „.] t (44) 

where the superscript T signifies the transpose and N the 
number of antenna elements. In this equation 5 f , is assigned 
to be unity. The square of the norm of Eq. (44) is 

|s| = I 8? 

i=0 

where M is the number of sources. The minimum norm method 
minimizes- this quantity and hence the name. The procedure is 
as follows: 
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1. Find the eigenvectors U of the covariance matrix R . 
Large eigenvalues correspond to signals and small 
eigenvalues correspond to noise. 

2. The vector D can be found from the signal subspace U s 
in Eq. (26) . 

3. The vector L 

subspace. The signal subspace U s can be partitioned as 



U s = 



can 


also 


be 


. subspace 


U, ' 


u oo 


u oi ■ 


U 0N-1 


u 10 


11 li 


U lN-l 


U M0 


U M1 ■ 


11 MN-1 






U c 



(45) 



where r\ H is the first row and U' s consists of the remaining 



rows of U s . These two matrices can be written as 



W = Wo 



U s = 



01 



The D vector can be obtained as 

1 



D = 



~Us r 1 / 



'V-r\ H r\) 



U 0N - /] 



U I0 U,j . . U 1N-I 



U M0 u Ml ■ ■ ' u MN-1 



(46) 



(47) 



4. Once the vector D is obtained, the angle estimator 
function can be written as 

1 



PmM = 



H 



a(e)DD H a(Q) 



( 48 : 
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So far several high resolution techniques have been 
applied to the problem of resolving the directions of 
arrival of incoming signals. From a statistical point of 
view, the signals arriving at the array can be regarded as 
random, and thus the signals and noise are uncorrelated. In 
the next chapter, we implement these methods based on finite 
observations . 
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IV. IMPLEMENTATION AND COMPUTER CODES 



A. PARAMETRIC DATA MODEL AND FLOWCHART 

The AOA estimation methods described in Chapter III 
were implemented in a computer code. Output data were 
gathered for a variety of arrays and signals, and a 
performance comparison made. This chapter describes the 
structure and operation of the computer simulation, which is 
written in Mat lab. The approach is model -based in the sense 
that it relies on certain assumptions made on the 
observation data. Since the natural causes responsible for 
signals and noise are often unrelated, it is customary to 
assume that the signals and noise are uncorrelated with each 
other. Figure 6 shows the flowchart of the computer program, 
with the main steps described in the following sections. 

1 . Input Parameters 

For an antenna array with N elements, the size of 
correlation matrix (correlation between i th and y'th 
element) is N by N . In this case, the maximum number of 
sources that can be detected cannot exceed (N- 1) . Because 
the effectiveness of these spatial spectrum methods are a 
function of signal-to-noise (SNR) ratio, we evaluate the 
performances of each antenna array by changing SNR at the 
input. Also, we must consider the spatial resolution 
capability of each estimation method when applied to each 
array geometry. For example, test signals that have the same 
azimuth or elevation angle are generally harder to resolve 
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Figure 6 . Computer simulation flowchart 

than those not located in a common principal plane. We also 
examine the effects of the number sampling points. The more 
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sampling points, the more accurate the AOA estimate but it 
occurs at the expense of increased processing time. 
Generally the antenna spacing must be less than 0.5 A in 
order to avoid the ambiguity problem. 

2 . Signal and Noise Models 



In practice, the white noise and signals are modeled as 
a independent and identically distributed random processes. 
That means we interpret the received signals as the sample 
of some waveform at certain specified instants of time, and 
is referred to as a realization of a discrete time random 
process. In all cases the signals that arrive can be 
regarded as random, we assume that signals s (n) are 

Laplacian probability distributed processes. Because all 
arrays are not diverse in polarization (i.e., they have all 
elements identically polarized) , we assume that the 

frequency and polarization of all received signals are the 
same. Multi-path and re-radiators in the vicinity of 
receiving site which may have a mutual coherence are not 
considered . 

The signal and white noise generation process are 

described by the following Matlab pseudo-code: 

u(n ) = rand(«)-0.5 

s(/?) = -(«(«)) * l°g(l - 2* | u{n) |) * A 



N„ 

M n) = -j=v( n) 



v(n) = randn(«) + j * randn( n) 



(49) 
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where n -0 , 1 , . . . , 4095 , A is the signal amplitude, N the 

peak amplitude of white noise and the noise v(n) is a 
complex Gaussian random variable with both the real and 
imaginary parts having a variance equal to 1. The Matlab 
functions rand and randn can be used to generate Gaussian 
distributed random values. The voltage ratio A / N 

determines the SNR; in dB it is 20 log l0 (A/N p ). For 
convenience we let A =1 so that 

N p =lO~ sm/20 . (50) 



3 . Antenna Array Model 



Based on the discussion in Chapter II the steering 
matrices for five array geometries can be determined. The 
steering matrix for the uniform linear array takes the form 



1 



1 



m = 



(51) 



e\p(j2n(N - l)d ■ • • e\p(j2n(N - l)d sin§ M ) 

where N is the number of antenna elements, M is the 
number of sources and is the AOA of the M th source, 




The steering matrix for a the non-uniform linear array 
is 
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A(^) = 



(52) 



_exp(j2nd N _, sin 4* 7 ) • • • exp( J /'2nfl r A ;. / sin<)) A/ ) 

where d t is the spacing between the reference element and 

the (/+l)th antenna. For a minimum redundancy array, the d l 
are specified in Table 1. 

The steering matrix for circular array is 
exp(j2nR a cos(-<j)/)sin9y) • ■ • Qxp(j2nR 0 cos(-^,)smQ M ) 






(53) 



exp(/2 nR 0 cos(4) ^ ) sin 9 , ) • • • exp(j2nR 0 cos(f VA/ ) sin 0 w ) 



where <j>„ OT = n§ 0 - <(>,„ , 9 m is the elevation angle of m th source, 

<j) m is the azimuth angle of m th source and R 0 is the 

radius of the circle. For the present analysis it is 

assumed that R a is equal to 0.52, 0<<j) m <2 n, and 0<9„, <^. 
The steering matrix for the spiral array is 

exp(j2nRj cos(-<j) / )sin9 / ) • • • zxp(j2nR, cos(-<j)/)sin9 A/ ) 



^4(9, (f>) = 



(54) 



Qxp(j2nR N cos(4> av ) sin 9 y ) • • • exp(j2nR N cos(<j) ) sin 9 A/ )_ 

where R n is the distance between n th antenna element and 
center point. In our model, the maximum R n is 0.52. 
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The steering matrix for the combination linear and 
circular array can be partitioned into sub-matrices 
associated with the linear and circular parts of the array 



2(0,<j>) 



1 

exp(y'27ic/sin0 / ) 
exp(y^7ic/sin0 / ) 
cxp(j2nd cos (-<(>,) sin 0, ) 



1 

exp(y'27tc/sin0 A/ ) 
exp(j4nds\r\Q M ) 
exp(J2nd cos(-<f> w )sin Q M ) 



exp(j2nd cos(<|> v/ ) sin 0 ; ) .... exp(y2^cos(<j) A , A/ )sin0 jW ) 

where d is 0.42. The first three rows in Eq. (55) are for 
the linear array elements. For these terms there is no 
dependence on arrival azimuth angle. 



4. Correlation Matrix and Spatial Spectrum Methods 

The output of the array can be calculated by using the 
signal, noise, and antenna models in Eq. (24) . First the 
Matlab function is used to find the singular value 
decomposition. The result is substituted into the various 
algorithms (Eq. (27), (41), (43), or (48)) with the 

appropriate antenna array steering matrix (Eq. (51), (52), 

(53) , (54) , or (55) ) . The outputs for MUSIC, PHD, minimum 

norm and Capon's beamformer are angle of arrival .vs. power. 
The power must be computed at every angle. We must steer the 
antenna array in all azimuth directions (0 to 360 degrees) 
and elevation direction (0 to 90 degrees) in order to find 
the AOAs . The output of ESPRIT is the direction of arrival 
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which can be obtained without a search. This method has only 
been applied to a uniform linear array. 

B. THE IMPLEMENTATION OF ARRAY MATRIX DECOMPOSITION 

The signal measured at the output of any element in the 
array differs from the signal actually received by an amount 
attributed to noise. Thus, we have an observed signal, which 
for each element of the array, consists of the actual 
received signal plus narrow-band noise. Calculations used 
£=4096 time sample points for each antenna element, so the 
size of received signal matrix is N by 4096, where N is 
the number of elements. The size of signal-in-space matrix 
is M by 4 096, where M is the number of sources. The 
observed signal matrix is equal to the received signal 
matrix times the propagation delay matrix. The observed 
signal matrix is the antenna output matrix G , from which we 
find the covariance matrix of the observed signal matrix 
defined by 

R = E [G * G T ] . 

The covariance matrix R has two important special 
properties. First, it is Hermitian, that is, it is equal to 
the conjugate transpose of itself. Second, it is positive 
semi-definite if element noise is present. Thus, for 
computer coding purposes 

R = E[G' G T ] = E [G G H ] . (56) 
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The size of R matrix is TV by N . The Matlab function SVD 
is used for singular value decomposition. It has the form 

[U, S, U] = SVD(R) 

where R = USU 11 (Eq. (26)) . The eigenvector U is an N by N 
orthogonal signal and noise subspace matrix, and S is an 
N by N diagonal matrix. The singular value decomposition 
is accurately computed by this function. The columns of U 

T 

are the eigenvectors of RR , and the entries along the 
diagonal of S are the correspondingly ordered nonnegative 
square roots of .the eigenvalues of R R , which are also the 

eigenvalues of RR . The numbers in S are the singular 
values of the matrix R , which is determined by the number 
of sources. After calculation of U and S', they are 
substituted into the spectral estimators in Eqs . (27), 

(43), and (48) to determine the AOAs . 

The MUSIC and PHD methods make use of the noise 
subspace matrix; the minimum norm method uses the signal 
subspace matrix. Capon's beamformer does not use the 
singular value decomposition of the covariance matrix. For 
2-D direction finding, it takes a significant amount of CPU 
time to search space in order to find the AOAs. The default 
search range in azimuth angle is 90 degrees. 
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C. GRAPHICAL USER INTERFACE 

To compare the relative performance of the various 
arrays and algorithms, the simulation must be run for many 
cases. A graphical user interface (GUI) was developed to 
vary the input parameters, as shown in Figure 7. Also 
included is an on-line help feature. The help button 
provides brief explanations of the input parameters. The GUI 
provides default values to demonstrate the angular 
resolution of MUSIC for a uniform linear array and a minimum 
redundancy array with different signal -to-noise ratios. 



antenna number: 
source number: 

Signal to Noise Ratio(dB): 

sample points: 



7 spacing: |. M 



CLOCK 



4. 



I 25 

I 4096 



CLOSE HELP 



1-D DIRECTION FINDING -90 ...90 

angle:| -is 1 -5 |, 5/1 is • 



2-D DIRECTION FINDING: EL 
angle:0...90 deg. AZ angle: 0... 



AZ angle: | 16 | 22 |; 28 : |. 34 
EL angle: j 10 | 15 -| 20 | 25 . 



2-D DIRECTION FINDING CIRCULAR ARRAY 



2-D DIRECTION FINDING SPIRAL ARRAY 



2-D DF LINEAR+CIRCULARARRAY 



LA-spacing: 0.4 
LA-ant No: 3 



Umfiorm and Non-uniform Linear Array 




MUSIC-resolution- MOVIE 


Linear Array with ESPRIT 



Antenna Geometry 

Uniform Linear Array 



Non-uniform linear array 



2-D DF Array 



Antenna Pattern 

Li near* circular array 



Circular array 



Figure 7 . Graphic user interface 



In the next chapter, we use this GUI to present a 
performance comparison of superresolution spectral 
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algorithms 
in several 



for elevating and analyzing their effectiveness 
antenna and emitter geometries. 
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V. 



SIMULATION RESULTS AND ANALYSIS 



In this chapter the results of the computer simulation 
are presented. The performances of the various algorithms 
with different array and emitter geometries are examined. 
Computer simulations were run for the conditions of: (1) two 
closely spaced AOAs for testing azimuth and elevation 
angular resolution, (2) various levels of additive white 
noise, (3) multiple signal sources, and (4) changes in the 
antenna element spacing. 

By comparing the test results, we can determine which 
combination of array and algorithm operates best under the 
given set of conditions. The power spectra are plotted as 
surfaces. For all cases, the z -axis represents magnitude of 
the power, the x-axis represents as azimuth angle and the 
y-axis the elevation angle. The majority of the cases use 7 
antenna elements; the exceptions are Figures 35 and 37 
which are for 10 elements. 

A. ANGULAR RESOLUTION OF ONE -DIMENSIONAL DF 

Initial simulations were run to determine the 
comparative performance of algorithms in resolving closely 
spaced AOAs. The signals were assumed to be of equal 
frequency, polarization and mutually incoherent. 

Simulations were run for two signals with angles of 
arrival of 15 and 18 degrees; the corresponding spectral 
estimates are shown in Figure 8. The rows correspond to 
spectra estimates for the PHD, MUSIC, minimum norm and 
Capon's beamformer algorithms, and the columns correspond to 
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the uniform array (left) and non-uniform linear array 
(right) . From this figure, we can see that the non-uniform 
array can achieve better angular resolution. All algorithms 
except the Capon's beamformer provide sharp distinguishable 
peaks, but MUSIC has a smoother spectrum (the separation of 
the peaks is not as pronounced) . Figure 9 shows an 
interference simulation for 6 sources. The non-uniform 
(minimum redundancy) array provides better dynamic range for 
all algorithms. 

B. TWO-DIMENSIONAL DIRECTION FINDING 

1. Azimuth Angular Resolution 

A range of SNRs were simulated. Figures 10, 11, and 12 
show the results for a high SNR (=25 dB) , at an elevation 
angle of 15 degrees and AAZ=5 degrees. The MUSIC, PHD and 
min-norm methods have similar results for all three 2-D 
array geometries except the Capon's beamformer. However, 
when the SNR is lower (15 dB) , the circular array performs 
better, as indicated by a comparison of Figures 13, 14, and 
15. Even though two peaks exist, there is no notch between 
them and therefore it may be difficult to correctly resolve 
the two spectral peaks. Figures 16, 17, and 18 show the 
results for SNR=30 dB, and two sources at the same elevation 
angle but separated in azimuth ( AAZ=3 degrees) . All arrays 
can resolve the two peaks. However, for the lower SNR level 
case, the circular array is generally superior. 
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Uniform linear array Non-uniform linear array 
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Capon MLM 
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AOA(degree) 



AOA(degree) 



Figure 8. Comparison of uniform (left) and non-uniform (right) 
linear arrays, 2 sources, 7 elements, separation angle 3 degrees 
f 15° and 18°) . SNR=25 dB 
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Figure 9. Unifrom (left) and non-uniform (right) linear arrays, 
SNR=25 dB, 7 antennas, 6 sources: 15°, 25°, 35°, 45°, 55°, and 65° 
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Figure 10. Linear+circular array, SNR=25 dB, 2 sources at 
(16°. 15°) and (21°. 15°). 7 elements 

MUSIC PHD 




EL(deg) 0 0 AZ(deg) 




EL(deg) 0 0 AZ(deg) 



Mini-norm 




EL(deg) 0 0 AZ(deg) 



Capon MLM 




50 
AZ(deg) 



Figure 11. Circular array, SNR=25 dB, 2 sources at (16°, 15°) 
and (21°. 15°) . 7 elements 
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Figure 12. Spiral array, SNR=25 
(21°, 15°), 7 elements 
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Figure 13. Circular array, SNR=15 dB, 2 sources at (16°, 15°) 
and (21°, 15°), 7 elements 
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Figure 14. Linear+circular array, SNR=15 cLB, 2 sources at 
(16°, 15°) and (21°, 16°), 7 elements 
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Figure 15. Spiral array, SNR=15 dB, 2 sources at (16°, 15°) 
and (21°, 15°) , 7 elements 
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Figure 16. Circular array, SNR=30 dB, 2 sources at (16°, 15°) 
(19°, 15°), 7 elements 
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Figure 17. Linear+circular array, SNR=30 dB, 
(16°. 15°) and (19°. 15°) , 7 elements 
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Figure 18. Spiral array, SNR=30, 2 sources at (16°, 15°) and 
(19°, 15°) , 7 elements 
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Figure 19. Circular array, SNR=25, 2 sources at (16°, 15°) and 
(16°, 18°), 7 elements 
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2 . Elevation Angular Resolution 



Simulations were run for the two signal case with 
angles of arrival at 18 and 15 degrees in elevation and 16 
degrees in azimuth. For the SNR equal to 25 dB, the 
corresponding spectral estimates are shown in Figures 19, 
20, and 21. The results for the circular and linear+circular 
arrays are almost the same, but the spiral array cannot 
correctly resolves the peaks. For a 15 dB SNR the data, 
which is shown in Figures 22 and 23, indicate that all 
algorithms, except Capon's beamformer, correctly resolve the 
two spectral peaks with a circular array. Only two 
algorithms can resolve the peaks for the circular array, so 
the elevation angular resolution of the linear+circular 
array appears to be better for this scenario. 

3 . Multiple Signal Sources 

Multiple signal simulations are shown in Figures 24, 
25, and 26. There are four signals arriving at azimuth and 
elevation angles of (10°, 20°) , (55°, 15°) , (30°, 16°) , and (45°, 25°) 
with equal amplitudes. The circular antenna array and 
linear+circular array can resolve the four peaks but the 
circular array has more sidelobes. Only the MUSIC algorithm 
correctly resolves the signal sources for the spiral array. 
In this scenario the linear+circular array is best. 
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Figure 20. Linear+circular array, SNR=25 dB, 2 sources 
(16°, 15°) and (16°, 18°), 7 elements 
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Figure 21. Spiral array, SNR=25 dB, 2 sources at (16°, 15°) 
and (16°, 18°), 7 elements 
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Figure 22. Circular array, SNR=15 dB, 2 sources at (16°, 15°) 
and (16°, 18°), 7 elements 



MUSIC 




EL(deg) 

Mini-norm 




0 20 40 60 80 

EL(deg) 



PHD 



60 




-20 1 * ' * - 

0 20 40 60 80 




EL(deg) 



Figure 23. Linear+circular array, SNR=15 dB, 2 sources at 
(16°, 15°) and (16°, 18°), 7 elements 
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Figure 24. Linear+circular array, SNR=25 dB, 7 elements, 4 
sources at (10°, 20°) , (55°, 15°) , (30°, 60°) , and (45°, 25°) 
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Figure 25. Circular array, SNR=25 dB, 7 elements, 4 sources at 
(10°, 20°) , (55°, 15°) , (30°, 60°) , and (45°, 25°) 
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Figure 26. Spiral circular array, SNR=25 dB, 7 elements, 4 
sources at (10° , 20°) , (55° , 15°) , (30°, 60°) , and (45°, 25°) 
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Figure 27. Circular array, SNR=15 dB, 7 elements, 2 sources 
at (16°, 15°), and (16°, 20°) 
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4 . Comparison of Azimuth and Elevation Angular 
Resolution 

The estimation methods and arrays were compared for 
two emitters with a 5 degree separation. Figures 27, 28, and 

29 are for SNR=15 dB and AEL=5 degrees; Figures 13, 14, and 

15 for SNR=25 dB) AAZ=5 degrees. We see that the elevation 
angular resolution is better than the azimuth angular 
resolution since the sidelobes in azimuth are higher than 
the sidelobes in elevation for all three array 

configurations. Even at lower SNRs, all three arrays 

employing the MUSIC method can distinguish the two sources 
in elevation. 

5. Improvement of Angular Resolution of Two-dimensional 
DF 

Next we investigate changing the radius of the circular 
array to 1.5 wavelengths. The results are shown in Figures 

30 and 31. Comparing Figures 13 and 14 with Figures 30 and 

31 (SNR=15 dB, AAZ=5 degrees) , the angular resolution is 

better with the larger radius as expected. Larger spacing 
may be susceptible to ambiguities. From Figures 32 and 33, 
where the range of azimuth angles is from 0 to 360 degrees, 
it is evident that there is no serious ambiguity problem. 
Thus, we can improve the angular resolution by . changing the 
antenna spacing. 
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Figure 28. Linear+circular array, SNR=15 dB, 7 elements, 2 
sources at (16°, 15°), and (16°, 20°) 
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Figure 29. Spiral array, SNR=15 dB, 7 elements, 2 sources 
at (16°, 15°), and (16°, 20°) 
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The circular array has a smoother spectrum and less 
spikes than linear+circular array, since the spacing of 

linear array part is larger than o . 5 X . From these 
simulation results, we conclude that the larger antenna 
spacing will enhance angular resolution at the expense of 
more ripples on the spectrum. We can suppress the sidelobes 
and enhance resolution by increasing the number of elements 
in the x-y plane and increasing the radius of circle and 



then adding a linear 


part 


on z 


-axis with a 


spacing less 


than 0.5 A, 


• 












Two 


elevation 


cases 


are 


considered . 


For 


a large 


elevation 


angle (80 


degrees) , 


there are 


four 


signals 


arriving 


at azimuth 


and 


elevation angles 


of 


(10°, 80°) , 



(90°, 80°) , (180°, 80°) , (270°, 80°) , with a SNR=15 dB . The 

simulations are shown in Figure 34 for the circular array (7 
elements, R n = 1.5 A), and in Figure 3 5 for the 

circular+linear array (3 elements, d= 0 . 4 A on z -axis) . In 
both cases the sidelobes increase in the elevation 
direction, but Figure 35 shows lower sidelobes. From these 
results we conclude that the elevation angular resolution is 
worse when the elevation angle is large. This is a well 
known result that is associated with the increase in 
beamwidth that occurs with array scanning [Ref. 9] . 

For small elevation angles (5 degrees) , there are four 
sources arriving at azimuth and elevation angles of 
(10°, 5°) , (90°, 5°) , (180°, 5°) , (270°, 5°) , and a SNR=15 dB . The 
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simulations are shown in Figure 36 for the circular array (7 
elements, R n = 1.5/1), and in Figure 37 for the 

circular + linear array (3 elements, cl = 0.4/1). Figure 37 
shows lower sidelobes in the plane of the sources. These 
sidelobes expand in the azimuth direction, implying that the 
azimuth angular resolution is worse when the elevation angle 
is small . This case is close to a broadside condition for 
the linear array and, consequently, the resolution ability 
is essentially that of the linear array. Thus the resolving 
ability has been improved at low elevation by adding a 
linear array part (3 elements on the z axis) . 

Based on the results of the computer simulations and 
performance comparisons, we present some conclusions and 
recommendations in the next chapter. 
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Figure 30. Circular array, R a = 1.5 Z, SNR=15 dB, 7 elements, 2 
sources at (16°, 15°), and (21°, 15°) 
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Figure 31. linear+circular array, d=1.5A., SNR=15 dB, 7 elements, 
2 sources at (16°, 15°), and (21°, 15°) 
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Figure 32. Circular array, R 0 =1.5A, SNR=15 dB, 7 elements, 2 
sources at (16°, 15°), and (21°, 15°), MUSIC method 
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Figure 33. Linear+circular array, d= 1.5 A, SNR=15 dB, 7 elements, 
2 sources at (16°, 15°), and (21°, 15°), MUSIC method 
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Figure 34. Circular array, SNR=15 dB, 7 elements, 4 sources at 
(10°, 80°) , (90°, 80°) , (180°, 80°) , and (270°, 80°) , MUSIC method 




Figure 35. Linear+circular array, SNR=15 dB, 10 elements , 4 sources 
at (10°, 80°) , (90°, 80°) , (180°, 80°) , and (270°, 80°), MUSIC method 
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Figure 36. Circular array, SNR=15 dB, 7 elements, 4 sources at 
(10°, 5°) , (90°, 5°) , (180°, 5°) , and (270°, 5°), MUSIC method 



40 - • 




EL(deg) u u AZ(deg) 

Figure 37. Linear+circular array, SNR=15 dB, 10 elements, 4 
sources at(10°,5°) , (90°, 5°) , (180°, 5°) , (270°, 5°) , MUSIC method 
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VI . CONCLUSIONS 



Over the past few decades superresolution techniques 
have been developed to resolve multicomponent plane wave 
fields. The most noteworthy of these is the MUSIC 
algorithm. Others include PHD, Minimum norm, ESPRIT and 
conventional and Capon's beamforming. Each of the methods 
has its advantages and disadvantages with regard to hardware 
and computational aspects (array configurations and element 
spacing) and emitter parameters (number of emitters and 
their spatial distribution) . 

This thesis has examined the performance of four of the 
estimation methods listed above for several linear, planar 
and volumetric arrays. Based on the formulations described 
in Chapter III, Matlab simulations of the methods were used 
to collect data for linear, circular, spiral and circular 
plus linear arrays. The emitter locations and spacings were 
varied, as was the signal -to-noise ratio. The data are 
presented as mesh surfaces, and the surface characteristics 
can be interpreted in terms of emitter resolution. 

For one -dimensional DF, the MRA provides the best 
dynamic range and resolution for a given number of elements 
and baseline. Since ESPRIT is generally applied to the 
spacing which is a known constant value, it is not directly 
applicable to MRAs . Using ESPRIT, the AOA can be obtained 
without a search technique, and in this regard is different 
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from the other four methods. Therefore, only the other four 
algorithms have been compared. 

The simulation results are in agreement with related 
two-dimensional DF performance simulations reported by 
Johnson and Miner [Ref. 10] for the MUSIC and Capon's MLM 
algorithms. The data in Chapter V leads to the following 
conclusions with regard to the circular and circular+linear 
array : 

(1) for high SNR levels, all four methods have about the 
same the performance, 

(2) elevation angular resolution is better than azimuth 
angular resolution for these three array configuration 
considered, 

(3) the azimuth angular resolution of the circular array is 
better than the spiral or circular+linear, 

(4) the elevation angular resolution of circular+linear 
array is better than spiral and circular+linear, 

(5) for multiple signal sources, the linear+circular array 
is better when antenna spacing is less than 0.5 wavelength, 

(6) larger element spacings result in better resolution, but 
there are more ripples and spikes; ambiguities will occur if 
the spacing becomes too large, 

(7) when the elevation angles of two signal sources are 
small, the azimuth angular resolution is worse than it is at 
large elevation angles, (because of slow falloff of the 
pattern in azimuth) 
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(8) when the elevation angles of two signal sources are 
large, the elevation angular resolution is worse than when 
elevation angles are small (because of slow falloff of the 
pattern in elevation) . 

From the data it appears that the MUSIC algorithm is 
more stable, almost always correctly resolves individual 
emitters, and more importantly, has low sidelobes. However, 
it has been demonstrated that under some conditions MUSIC 
fails to provide the correct number of peaks (Ref. 11] . 

The circular array has best azimuth resolution but the 
results of linear+circular array are close to that of the 
circular array. Based on multiple signal sources the 
linear+cicular array is probably the better choice. When the 

antenna spacing is larger (say R 0 = 1.5 wavelength), the 

circular array is better. Having more elements on a larger 
diameter does provide some advantages. The distance between 
two adjacent antennas will be small to avoid high sidelobes 
and offers better azimuth resolution, especially for small 
elevation angle. 

Since we use subspace methods to figure out the AOA, 
the number of sources must be less than the number of array 
elements. Future work should consider how to deal with this 
limitation as well as improvements to the low signal-to- 
noise-ratio case. Any techniques that would allow a 
decrease in the search time for 2-D direction finding would 
also be valuable. 
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APPENDIX A 



THESIS MAIN PROGRAM 



clear all 
close 

cv=version; 
if str2num (cv ( 1 ) j <5 
figure ( 1 ) , 
uicontrol ( ... 

’ Style 1 , 1 text ' , ... 

' Units ' , ' normalized ' , ... 

' Position ' , [0.01 0.3 0.9 0.1 ], ... 

' BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ' , ' left ' , ... 

’String', 'This program cannot be run on MATLAB 4. X, please use 
MATLAB 5 . X ’ ) ; 

pause (5), close,exit 



end 

load nsource 
msens=7; 
snr=25 ; 

%nsource=2 ; 

uicontrol (' Units ',' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' FontSize ' , 10 , ... 

' ForegroundColor ' , [ 0 01], ... 

'Position', [23.5 288.75 125 17.25], ... 

’ String ' , ' antenna number : ' , ... 

’ Style ’ , ' text ' , ... 

'Tag' , ' StaticText 1 ' ) ; 
uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

’ FontSize ' , 10 , ... 

' ForegroundColor ',[ 1 00], ... 

' Position’ , [22 . 75 273.75 125.75 15], ... 

' String ' , ' source number : ' , ... 

'Style' , 'text' , ... 

' Tag ' , ' StaticTextl ' } ; 
uicontrol (' Units ',' points ' , ... 

' BackgroundColor ', [0 . 8 0.8 0.8], ... 

’ FontSize ' , 9, ... 

’ Position' , [22 254 125 19], ... 

' String ’,' Signal to Noise Ratio(dB):', ... 
' Style ' , ' text ' , ... 

'Tag' , ' StaticText2 ' ) ; 



msensl = uicontrol (' Units points ' , ... 

' BackgroundColor ',[ 0 10], ... 

'Callback' , ' inn3 ' , ... 

' FontSize ' , 12, ... 

'Position' , [150.5 290.25 27 15], ... 

' String ', num2str (msens ) , ... 
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’Style’ , ’edit 1 , . . . 

’Tag' , ' EditText 1 ' ) ; 

nsourcel = uicontrol (' Units points ’ , ... 

’ BackgroundColor ' , [ 0 10], ... 

' Callback resl ' , ... 

' FontSize’ ,12, ... 

' Position’ , [150. 5 274.5 27 15], ... 

' String ’, num2str (nsource) , ... 

' Style ' , ’edit ’ , ... 

' Tag ' , ' EditText 1 ' ) ; 

snrl = uicontrol (' Units points ’ , ... 

' BackgroundColor ’, [0 10], ... 

’Callback’ , ’ inn3 ’ , . . . 

' FontSize ' , 12, ... 

’Position’ , [150.5 258 27 15], ... 

’String ’ , num2str (snr ) , ... 

' Style ' , ' edit ? ) ; %, ... 

uicontrol {’ Units points ’ , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' FontSize ' , 10, ... 

' ForegroundColor ' , [ 1 00], ... 

’ Position' , [40. 9655 242.69 100.552 17.3793], ... 

' String ',’ sample points:’, ... 

’Style' , ’text’ ) ; 

bnsamp=uicontrol { ’ Units ' , ' points ' , ... 

' BackgroundColor ’,[ 0 10], ... 

' Callback inn3 ' , ... 

’ Position’ , [150.207 242.069 26.069 14.2759], ... 

’String’ , ’4096’ , ... 

’Style' , ’edit’ ) ; 

uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' FontSize ’ , 10, ... 

' ForegroundColor ', [0 01], ... 

’ Position’ , [178.759 287.379 54 16.7586], ... 

’ String 1 , ’ spacing : ' , ... 

’Style’, ’text'); 

bdspace= uicontrol (' Units ’,’ points ’ , ... 

’ BackgroundColor ’,[ 0 10], ... 

’Callback’ , ’ inn3 ’ , ... 

’Position’, [232.138 288.621 26.069 14.2759], ... 

’String’ , '0.4 ’ , ... 

' Style ' , ' edit ' ) ; 
aw=clock; 

tname = [’Direction Finding Analysis year:’, num2str (aw { 1 ) ) , ' , 

month: ', num2str (aw ( 2 ) ) ,’, — day: ’, num2str (aw { 3) ) ] ; 

set(gcf, ' name ’, tname) 

nsource=str2num ( get (nsourcel , ' string ’ ) ) ; 
bdspacel= uicontrol (' Units ’, ’points ' , ... 

’ BackgroundColor ',[ 0 10], ... 

' Callback ’,' inn33 ' , ... 

' Position’ , [252 27 27 12], ... 

’String' , ’0.4 ' , ... 

' Style ’ , ' edit ' ) ; 

bdl= uicontrol (’ Units ’,' points ' , ... 

' BackgroundColor ',[ 0 10], ... 



66 



' Callback’ , ’ inn 3 3 ’ , ... 

’ Position' , [252 14 27 12], ... 

’String’ , ’ 3 ' , ... 

' Style ’ , ' edit ' ) ; 



uicontrol (’ Units points ' , ... 

' BackgroundColor ' , [0 . 8 0.8 0.8], ... 

' FontSize ’ , 10, ... 

’ ForegroundColor ' , [ 0 01], ... 

’ FontWeight ' , ' bold ' , ... 

’Position’ , [11. 1724 222.828 218.25 15.7172], ... 

' String ’ , ' 1-D DIRECTION FINDING -90 ...90 deg’, ... 
' Style ’ , ’ text ’ , ... 

' Tag ’ , 1 StaticText3 ' ) ; 

uicontrol (' Units ’,’ points ’ , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' FontSize ',10, ... 

' ForegroundColor ' , ' r ' , ... 

’ FontWeight ' , 'bold' , ... 

’ Position’ , [16.25 207.5 41.75 15], ... 

’ String ' , ’ angle : ’ , ... 

' Style ' , ' text ’ , ... 

' Tag ' , ’ StaticText4 ' ) ; 



for n=l:8 

bbl=uicontrol ( ’ Units ' , ’ points ' , ... 

'BackgroundColor ', [0. 8 0.8 0.8], ... 

’ Position ’ , [56.25+(84-56.25) * (n-1) 207 26.25 14.25] 
’ String ' , ' ’ , ... 

’ Style ' , ' text ' ) ; 

end 

%clear bbl 
%bbl= [ ] 

for n=l:nsource 

bbl (n) =uicontrol ( ’ Units ’ , ’ points ' , ... 

' BackgroundColor ', [0 10], ... 

’ Callback ’ , ’ inn3 ' , ... 

’ Position' , [56.25+ (84-56.25) * (n-1) 207 26.25 14.25] 
’ String 1 , int2str (-25 + 10*n) , ... 

' Style ' , ’ edit ' ) ; 

^nd 

size (bbl ) ; 
bbl; 

uicontrol (' Unit s ',’ points ' , ... 

’ BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' FontSize ’ , 9 , ... 

’ ForegroundColor ’ , ' r ' , ... 

' FontWeight ’ , ’ bold ' , ... 

'Position' , [2.75 126.75 58.75 12.75], ... 

'String', 'AZ angle:', ... 

' Style ' , ' text ' ) ; 



for n=l:8 
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bb2(n) = uicontrol (' Units ',' points ' , ... 

' BackgroundColor ' , ( 0 . 8 0.8 0.8], ... 

' Position' , [58.25+ (90-62.25) Mn-1) 126.75 26.25 14.25], ... 

' String 1 , ' ' , ... 

' Style ' , ’ text ' ) ; 



end 

clear bb2 
bb2= [ ] ; 

for n=l : nsource 

bb2(n) = uicontrol (' Units ', 1 points ' , ... 

' BackgroundColor ', [0 10], ... 

' Callback ' , ' inn3 ' , ... 

’Position', [58.25+(90-62.25) *(n-l) 126.75 26.25 14.25], ... 
1 String ’, num2str ( 10+6*n) , ... 

' Style ' , 1 edit ' ) ; 



end 

uicontrol (' Units ',' points ' , ... 

' BackgroundColor ', [0 . 8 0.8 0.8], ... 

’ FontSize ' , 10, ... 

' ForegroundColor 1 , [ 0 01], ... 

1 FontWeight ' , ’ normal ' , ... 

’ Position' , [12 156 207 32.25], ... 

•String', '2-D DIRECTION FINDING: EL angle:0...90 deg. AZ angle 
0. . . 360 deg’ , ... 

' Style ' , ' text ' ) ; 



uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' FontSize ' , 9, ... 

' FontWeight ' , ' bold ' , ... 

' Position' , [2 108 58.75 12.75], ... 

' String ' , ' EL angle : ' , ... 

' Style ' , ' text ' ) ; 
for n=l:8 

bb3(n) = uicontrol (' Units ',' points ' , ... 

' BackgroundColor ', [0 . 8 0.8 0.8], ... 

' Position’ , [58+ (90.75-63) * (n-1) 108 26.25 14.25], ... 

'String' , ' ' , ... 

'Style' , 'text ' ) ; 

end 

clear bb3 
bb3= [ ] ; 

for n-1 : nsource 

bb3(n) = uicontrol (' Units ',’ points ' , ... 

' BackgroundColor [0 10], ... 

' Callback ' , ' inn3 ' , ... 

' Position' , [58+ (90.75-63) * (n-1) 108 26.25 14.25], ... 

'String ' , num2str (5 + 5*n) , ... 

' Style ' , ' edit ' ) ; 

end 

bval= uicontrol (' Units ',' points ' , ... 

' BackgroundColor ', [0 0 1], ... 

' Callback' , ' inn3 ' , ... 
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' FontSrze ' , 10, ... 

f ForegroundColor ' , [ 1 10], ... 

1 Position' , [190. 552 153.31 36.6207 19.2414], ... 

’String' ,[' 90 ’ ; ’ 180 ' ; ’ 270 ’ ; ’ 360 ’ ] , ... 

' Style ' , ' popupmenu ' , ... 

' Tag ’ , ' PopupMenul ' ) ; 

clear bearingll 
clear bearingn 
clear ang 

msensll=str2num (get ( msensl , ’ string ' ) ) ; 
snr=str2num(get (snrl , * string ’ ) ) ; 
dspace=str2num (get (bdspace, ’ string ' ) ) ; 
nsamp=str2num ( get (bnsamp, ’ string ' ) ) ; 
thrv=get (bval , 'Value') 
if ( thrv==l ) 
thr=90 ; 

elseif (thrv==2 ) 
thr=l 8 0 ; 
elseif (thrv==3) 
thr=27 0 ; 
elseif ( thrv==4 ) 
thr=360 ; 

end 

for n=l:nsource 

bearingll (n) =str2num (get (bbl (n) , 1 string ’ ) ) ; 

bearingn (n) =str2num (get ( bb2 (n) , ' string ' ) ) ; 
ang (n) =st r2num ( get (bb3 (n) , ' string ’ ) ) ; 

end 

bearingll 

bearingn 

ang 

uicontrol (' Units points ' , ... 

'BackgroundColor ', [0.752941 0.752941 0.752941], 

' ButtonDownFcn’ , ’ ctlpanel SelectMoveResize ’ , ... 

'Position' , [288 151.5 125.25 69], ... 

'Style' , ' frame' , ... 

' Tag ' , ’ Frame 1 ' ) ; 

uicontrol (' Units ',' points ' , ... 

' ButtonDownFcn ',' ctlpanel SelectMoveResize', ... 

' Callback ',' doam ' , ... 

' FontSize ' , 9, ... 

’ FontWeight ' , ' bold ' , ... 

'Position' , [230.25 224.25 182.25 28.5], ... 

' String ',' Unifiorm and Non-uniform Linear Array' 

' Tag ' , ' Pushbuttonl ' ) ; 
uicontrol {' Units ',' points ' , ... 

' ButtonDownFcn ',' ctlpanel SelectMoveResize ' , ... 

' Callback ',' eprit ' , ... 

' FontSize ' , 9, ... 

' FontWeight ' , 'bold' , ... 

'Position' , [293.25 156 116.25 27], ... 

'String',' Linear Array with ESPRIT', ... 

' Tag ' , ' Pushbutton2 ' ) ; 
uicontrol (' Units ',' points ' , ... 

'ButtonDownFcn', 'ctlpanel SelectMoveResize', ... 
'Callback' , 'mvml ' , ... 
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' FontSize ' , 9, ... 

' Font Weight ' , ' bold ' , ... 

1 Position 1 , [292.5 186.75 117.75 28.5], ... 

’ String ’ , ’ MUSIC- re solution- MOVIE 1 , ... 

' Tag ’ , 1 Pushbutton3 1 ) ; 
uicontrol ( 1 Units ', 1 points ' , ... 

' ButtonDownFcn ’ , ' ctlpanel SelectMoveResize ’ , . 

1 Callback ',' doag4 ’ , ... 

' FontAngle * , ' italic ' , ... 

' FontName’ , ’TIMES’ , ... 

' FontSize ’ , 10, ... 

’ FontWeight ' , 1 bold ’ , ... 

' Position’ , [10.5 72.75 270.75 24.75], .. 

’String’ , ’2-D DIRECTION FINDING CIRCULAR ARRAY 
' Tag ’ , ’ Pushbutton4 ’ ) ; 
uicontrol ( ’ Units ',' points ' , ... 

' ButtonDownFcn ctlpanel SelectMoveResize', . 

’ Callback ' , ’ spira ' , ... 

' FontAngle ’ , ’ italic ' , ... 

' FontName ' , 'TIMES' , ... 

’ FontSize ' , 10, ... 

' FontWeight ' , ' bold ' , ... 

' Position' , [11.25 41.25 270 27], ... 

'String' , '2-D DIRECTION FINDING SPIRAL ARRAY’, 

' Tag ' , ' Pushbutton5 ' ) ; 
uicontrol ( ' Units point s ' , ... 

' ButtonDownFcn ’,' ctlpanel SelectMoveResize’, . 

. 'Callback' , ' cl ’ , ... 

' FontAngle ' , ' italic ' , ... 

' FontName ' , ' TIMES' , ... 

1 FontSize ' , 10, ... 

' FontWeight ' , ’bold' , ... 

'Position' , [10.5 14.25 190 22.5], ... 

’String' 2-D DF LINEAR+CIRCULAR ARRAY’, ... 

’ Tag ' , ' Pushbutton6 ' ) ; 
uicontrol (' Unit s points ' , ... 

' BackgroundColor ' , [0 . 8 0.8 0.8], ... 

’ FontSize ' , 9, ... 

’ ForegroundColor ' , [0 01], ... 

’Position', [200.759 22.25 54 16.7586], ... 

’ String ' , ' LA- spacing : ' , ... 

'Style' , 'text' ) ; 

uicontrol (' Units ',’ points ' , ... 

' BackgroundColor ',[0.8 0.8 0.8], ... 

' FontSize ' , 9, ... 

' ForegroundColor' ,[001], ... 

'Position', [200.759 8.25 54 16.7586], ... 

’ String ' , ’ LA-ant No : ' , ... 

'Style' , 'text' ) ; 

uicontrol ( ' Units ' , 'points ' , ... 

'ButtonDownFcn', 'ctlpanel SelectMoveResize', . 
'Callback ' , 'close ' , . . . 

' FontSize ' , 12, ... 

’ FontWeight ' , ' bold ' , ... 

'Position', [291.75 258 52.5 25.5], ... 

'String' , 'CLOSE' , ... 
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f Tag ' , f Pushbutton? ' ; ; 



uicontrol ( ' Units ' , ’points * , ... 

'ButtonDownFcn', * ctlpanel SelectMoveResize', 
’Callback’ , ’ ! clock’ , . . . 

’ FontSize ’ , 12 , ... 

' FontWeight ' , ' bold ’ , ... 

'Position*, [291.75 284.25 102.75 24.75], ... 

'String' , 'CLOCK' , ... 

' Tag ' , ' Pushbutton7 ' ) ; 
uicontrol ( ' Units points ' , ... 

' BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

' But tonDownFcn ' , 'ctlpanel SelectMoveResize ' , 

' FontName ' , ' times ' , ... 

'FontSize' , 13, ... 

' FontWeight ' , ' bold ' , ... 

' ForegroundColor ’ , [ 1 00], ... 

’ Position’ , [287 .25 127.5 128.25 23.25], 
' String ’,' Antenna Geometry', ... 

' Style' , ’ text ' , ... 

’ Tag ’ , ' StaticTextS ' ) ; 
uicontrol ( ' Unit s ',' point s ' , ... 

' ButtonDownFcn ' , ' ctlpanel SelectMoveResize ' , 

' Callback ',' sp2 ' , ... 

' FontSize ' , 9, ... 

' FontWeight ' , ’ bold ' , ... 

'Position' , [291 87.75 122.25 22.5], ... 

’ String' ,' Non-uniform linear array', ... 

'Tag' , ' Pushbutton8 ’ ) ; 
uicontrol ( ' Units ',' points ' , ... 

' ButtonDownFcn ' , ' ctlpanel SelectMoveResize ' , 
'Callback' , 'spl ’ , ... 

' FontSize ' , 9, ... 

' FontWeight ' , 'bold' , ... 

'Position' , [290.25 111.75 122.25 22.5], 
' String ',' Uniform Linear Array', ... 

' Tag ' , ' Pushbutton8 ' ) ; 
uicontrol ( ' Units ',' point s ' , ... 

'ButtonDownFcn', 'ctlpanel SelectMoveResize', 
'Callback' , *g5' , ... 

’ FontSize ' , 10 , ... 

' FontWeight ' , 'bold' , ... 

' Position’ , [291 64.5 122.25 22.5], ... 

’ String ' , ' 2-D DF Array ' , ... 

' Tag ' , ' Pushbutton8 ' ) ; 
uicontrol’ ( ' Units ',' points ' , ... 

'ButtonDownFcn', 'ctlpanel SelectMoveResize', 

' Callback ' , ' helpwin ( str , ttl ) ' , ... 

' FontSize ' , 15, ... 

' FontWeight ' , 'bold' , ... 

' Position ' , [ 34 6 . 5 258.75 47.25 24], ... 

' String ' , ' HELP ' , ... 

’ Tag ' , ' Pushbutton8 ' ) ; 
uicontrol (' Units point s ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

'ButtonDownFcn', 'ctlpanel SelectMoveResize', 

' FontName ' , ' times ' , ... 
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' FontSize ' , 13 , ... 

' FontWeight ' , ' bold ' , ... 

' ForegroundColor ' , [ 1 00], ... 

' Position' , [289. 5 41.25 128.25 23.25], ... 

' String Antenna Pattern', ... 

'Style', 'text' ) ; 
uxcontrol (' Units ',' points ' , ... 

' ButtonDownFcn ' , ' ctlpanel SelectMoveResize ' , ... 

' Callback ' , ' cclap ' , ... 

' FontSize ' , 9, ... 

' FontWeight ' , 'bold' , ... 

’Position' , [292.5 27.75 122.25 22.5], ... 

'String',' Linear+circular array'); 
uicontrol ( ' Units ' , 'points ' , ... 

' ButtonDownFcn ',' ctlpanel SelectMoveResize', ... 
’ Callback clap ' , ... 

' FontSize ' , 9 , ... 

' FontWeight ' , 'bold' , ... 

•Position’ , [292.5 3.75 122.25 22.5], ... 

' String ',' Circular array’); 



str 



This program was created by LIN, KU-TING, TAIWAN ' 
Date : 9/ 98 



geometry 



process 

multipl 

predict 



mg ' 



Thesis topic: ' 

Comparsion of super resolution algorithms with different array 

i 

for radio ‘direction finding' 

i 

Thesis adviser : Prof . David C.Jenn' 
second reader: Prof. Pace Phillip' 



We evaluate the ef fectiveness of various array designs and 



method for two-dimensional direction finding and resolving 

i 

emitters. Using five different spatial spectrum methods to 
the ' 

direction of arrival ( DOA) ' 



Discussion of Topics:' 

1. spatial spectrum estimation methods' 

(a) MUSIC (multiple signal classif ication ' 

(b) ESPRIT' 

(c) PHD ( Pisarenko harmonic decomposit ion) ' 

(d) minmum-norm* 

(e) Capon beamformer ' 

2. Array Geometry' 

(a) uniform linear array(l-D)' 

(b) non-uniform linear array(l-D) ' 

(c) circular array(2-D)' 

(d) spiral array(2-D)' 

(e) linear+circular array (2-D) ' 

3. Antenna Pattern' 

(a) circular array' 

(b) linear+circular array' 

4. Angular Resolution' 



72 



HOW TO USE?' 



array 1 



results ! ' 

» 

length) ' 



Input parameter' 

(a) antenna number: less than 10 for non-uniform linear 

(b) source number: less than antenna number* 

(c) SNR (dB) : signal to noise ratio.' 

(d) sample points: the more points, the more accurate 

( e ) spacing : 1 -D DF antenna spacing (normalized to wave- 



element 
axis ) of 
axis ) of 



menu ! 



2-D DF the distance from center point to 

(X-Y plane) ' 

( f ) LA spacing : antenna spacing for linear array part (Z- 
linear+circular array’ 

(g) LA-ant No:antenna number for linear array part(Z- 
linear+circular array’ 

(h) l-D direction finding angle range: 0^180 deg’ 

(i) 2-D DF EL angle:0-90 degree’ 

( j ) 2 — D DF AZ angle:0~360 degree' 

note .-choose the AZ searching range --> push popup 

only four choices : 90 , 1 80 , 27 0, 360 degree. 



moment ! 



2. When you change the ’’source number ”, please wait one 

3. Check antenna array geometry and antenna pattern' 

4. Choose 1-D or 2-D DF --> Push button, See what you get!’ 



’ * >--Attention--<* ’ 

’ 1. The 2-D DFs take some time to search the AOAs, ’ 

' please be patient to wait, the word "HOLD ON!!" will show 

up ! ' 

’ 2. ESPRIT method just for uniform linear array: 20 runs' 

' 3. MUSIC-resolution-movie DEMO — > four tagerts fixed, one 

moves . ' 

' sample points=200' 

’ 4. This program CANNOT be run on Mat lab 4.X, at least 



Matlab 5.X 

i 

10 (=7+3) ! 

» 

resolution 

i 

sidelobe ’ 

i 

f 

sidelobe ’ 



5. The default antenna number of linear+circular array is 

6. The larger radius of circular array, the better 

7. When the elevation angle of AOA is small (5 degree), the 

expand on the azimuth direction. ’ 

8. When the elevation abgle of AOA is large (80 degree), the 



' expand on the elevation direction.’ 

’ 9. When the radius of circular array is larger (1.5 

wavelength ) , there ’ 

' are more spi kes, supressing the spikes by adding linear 

array on the Z-axis’ 

' ’ } ; ttl = ’ DOA ' ; 

%helpwin (str, ttl) ; 

text (0.65, 1.02, set str ( 108 ) , ’ font name ' , ’ symbol ’ , ’ font size ’ , 20) 
axis ( ’ of f ’ ) 
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APPENDIX B 



SPECTRUM CALCULATION FOR CIRCULAR ARRAY 



format compact 
bearing=bearingn; 
msens=msensl 1 ; 
msource=nsource; 
def_msens = 8; 
def_msource = 2; 

smat = zeros ( nsamp, msource) ; for -n = l:msource 
u = rand (nsamp, 1 ) -0.5; 

smat(:,n)= -sign(u) .* log ( l-2 + abs (u) ) ; 

end 

< > 

bearing = bearing * pi / 180; % degrees to rad 

phi = 2 4 'pi/msens ; 
al= [ 0 : msens-1 ] '‘'phi; 
for nnn=l : msource 
cira (nnn, : ) =ql ; 

end 

for nb=l:msens 

cb ( : , nb) =bearing ’ ; 

end 

amat2= ( (2*pi*dspace*cos (cira-cb) ) ) ’ ; 
for ns=l:msource 

a 3 ( : , ns) =amat2 ( : , ns) *sin (ang (ns) * pi/ 180) ; 

end 

amat 3=exp ( j *a3 ) ; 

% the noisefree sensor signal 
smat 2=smat + amat 3 . ’ ; 

[nsamp, msens] = size(smat2); 
nsource=msource ; 

%snr=input ( ’ snr= ’ ) ; 

namp=l/ (sqrt (2) - r 10 A (snr/20) ) ; 

for n=l : msens 

nmat ( : , n) =randn (nsamp, 1 ) + j *randn (nsamp, 1 ) ; 
end 

Inmat = nmat - ones ( nsamp, 1 ) * mean (nmat); 

smat2=smat 2+namp*nmat ; 
ymat=smat2 ; 

% estimate the spatial cross-cumulant matrix 

cmat = zeros (msens , msens ) ; 

ymat = ymat - ones (nsamp, 1 ) * mean(ymat); % remove mean 

cmat - conj(ymat’ + ymat) / nsamp; % correlation matrix 

% determine number of sources 

[umat, smat, vmat] = svd(cmat); 
svec = diag(smat); 

format compact 



estimate bearing spectra 
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prmat = zeros (msens , msens ) ; 
for i=l:msens 

i f ( svec ( i ) > 0 ) 

prmat = prmat + vmat ( : , i ) + vmat ( : , i ) ' / svec ( i ) ; 

end 

end 

prvec = prmat(:,l); 
dtheta=l ; 

theta = [ 0 : dtheta : thr ] ' * (pi/180); 

va= [0:1: 90] * (pi/180) ; 
pp= (length (theta) ) ; 
for nnl=l:pp 

ciral (nnl, : ) =ql; 

end 



for nbl=l:msens 

cbl ( : , nbl) =theta; 



end 



omega = 2*pi*dspace*cos ( ciral-cbl ) ; %< 
@@@@@@@@*** **%%%% AAAA 



mth = length (theta) ; 

mlc = zeros (mth, 1 } ; 



gvec = umat ( 1 , 1 : nsource ).* ; 
gmat = umat (2 : msens, 1 : nsource) ; 

arvec = [1 ; - conj(gmat) * gvec / ( 1 - gvec’ *gvec) ]; 



tmat = vmat (:, nsource+1 : msens) ; 

Iwvmat = tmat * diag (ones (msens-nsource, 1 ) . / svec ( nsource + 1 : msens ) ) * 

tmat 1 ; 

tmat = tmat * tmat ’ ; 
nv=length ( va ) ; 
if thr==360 
for bb=l:nv 
for k =l:mth 

steer = exp(j * omega ( k, :)* *sin (va (bb) )) ; 
smusl(k,bb) = 1./ (steer* * tmat * steer); 
end 
end 



else 

for bb=l:nv 



for k =l:mth 



steer = exp ( j 
beaml ( k, bb) = 1 . 
sminl ( k, bb) = 1 . 
smusl ( k, bb) = 1 . 
spisl ( k, bb) = 1 . 



* omega ( k, :)* *sin (va (bb) )) ; 

/abs (steer* + inv( cmat)* steer) ; 

/ abs (arvec.* * steer) A 2; 

/ (steer* * tmat * steer); 

/ abs (steer* * vmat (:, msens ) ) A 2; 



end 

end 

end 
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val= [ -90 : 1 : 90] + (pi/ 180 ) ; 
mtl=length { val ) ; 

display the angular spectra 

theta - theta* 180/pi ; 
va=va*180/pi; 
val=val* 180/pi; 

top=0 . 95 ; spacing=0 . 02 ; 
left=0. 80; 
btnWid=0 . 1 ; 
btnHt=0 . 04 ; 
bottom=0 . 05 ; 

figure ( 1 ) 

uicontrol (' Units ',' points ' , ... 

1 BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

' FontSi ze ' , 1 9, ... 

' ForegroundColor 1 , [ 0 . 8 0.8 0.8], ... 

' Position ', [178 264 100 19], ... 

'String' , 'HOLD ON ! ! ' , ... 

' Style ' , ' text ' , ... 

' Tag ' , ' StaticText2 ' ) ; 
f a=get ( 0, ' current figure ' ) ; 
figure ( f a+1 ) , 

as=10 *logl0 (abs ( ( smusl+eps) ' ) ) ; 

asl=10*logl0 (abs ( (smusl+eps) ' ) ) ; 

%mesh ( theta ' , va , 1 0 * log 10 (abs ( ( smusl+eps )'))); 

bb=10*logl0 (abs ( (smusl+eps) ' ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

if thr==360 

mesh ( theta ', va, bias+10*logl0 ( abs ( (smusl+eps)'))); 

%title ( [ 'MUSIC, circular array, SNR= ' , num2str (snr, 3) , ' , radius-C- 

spacing : ' , num2str (dspace, 3) ] ) ; 

zlabel ( ' power (dB) ' ) ; 

xlabel ( ' AZ (deg) ' ) 

ylabel ( ' EL (deg) ' ) 

axis([0 90*thrv 0 90 -20 as] ) ; 

%xlabel ( [ 'ANTENNA NUMBER= ' , int2str (msens) ] ) ; 

%ylabel ( [ ' source number=' , int2str (nsource) ] ) ; 

zlabel ( ' power ’ ) 

break 

else 

as2=10*logl0 (abs ( (spisl+eps) ' ) ) ; 
as3=10~logl0 (abs ( (sminl+eps) ' ) ) ; 
as4=10 *logl0 (abs ( (beaml+eps) ' ) ) ; 

subplot ( 221 ), mesh ( theta',va,bias+10*logl0(abs( ( smusl+eps ) ' ) ) ) 
zlabel ( ' power (dB) ' ) ; 
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title ( 'MUSIC' ) ; xlabel ( ’ AZ (deg) ’ ) 
ylabel ( ’ EL (deg) ' ) 

axis([0 90 + thrv 0 90 -20 as]); 

-oxlabel ( [ ’ ANTENNA NUMBER^ ' , int2str (msens ) ] ) ; 
o ylabel ( [ ' source number= 1 , int2str (nsource) ] ) ; 
o zlabel ( ' power (dB) ’ ) 

bb=10*logl0 (abs ( (spisl + eps) ' ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

subplot ( 222 ) , mesh ( theta ',va,bias+10*logl0(abs( ( spisl+eps )'))); 

title ( ' PHD' ) ; 

zlabel ( ' power (dB) ’ ) ; 

xlabel ( ' AZ (deg) ' ) 

ylabel ( 'EL (deg) ' ) 

axis([0 90*thrv 0 90 -20 as]); 

%xlabel ( 'Az-angle' ) ; 

bb=10*logl0 (abs ( (sminl+eps) ' ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

subplot (223) ,mesh ( theta ',va, bias+10*logl0 (abs ( (sminl + eps) '))); 

title ( ' Mini -norm ' ) ; 

zlabel ( 'power (dB) ' ) ; 

xlabel ( 'AZ (deg) ' ) 

ylabel ( 'EL (deg) ' ) 

axis([0 90*thrv 0 90 -20 as] ) ; 

%ylabel { ' El-angle ' ) ; 
bb=10^1ogl0 (abs ( (beaml+eps ) 1 ) ) ; 
bias=min (min (bb) ) ; 

%as=max (max (bb) ) ; 

subplot (2 24 ) ,mesh( theta ',va,abs(bias)+10*logl0 (abs { (beaml + eps) ' ) ) ) 

title ( 'Capon MLM ' ) ; 

zlabel { 'power (dB) ' ) ; 

xlabel ( ' AZ (deg) ’ ) 

ylabel ( ' EL (deg) ' ) 

end 

AZ=bearing* 130 /pi ; 

EL=ang; 



%t it le ( [ 'ANTENNA NUM3ER= ’ , int 2s tr (msens ) , ' SOURCES 
NUMBER= ' , int2str (msource) ] ) 

% xlabel ( [ ' AZ= ' , int2str (AZ { 1 : msource) ) , ' degree ' ] ) 

% ylabel ( [ ' EL= ' , int2str (EL ( 1 : msource) ) , ' degree ' ] ) 
tname = [' CIRCULAR ANTENNA ARRAY, SNR= ' ,int2str ( snr )] ; 

set (gcf , ' name ' , tname) 
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APPENDIX C - SPECTRUM CALCULATION FOR SPIRAL ARRAY 



format compact 
bearing=bearingn; 
def_msens = 8; 
def_msource = 2; 
msource=nsource ; 
msens=msensll ; 

smat = zeros (nsamp, msource ) ; 
for n = 1: msource 
u = rand (nsamp, 1 ) -0.5; 

smat (:,n)= -sign(u) . * log { l-2*abs (u) ) ; end 
% the propagation matrix 
bearing = bearing * pi / 180; 

phi =2* pi /msens ; 

ql= [ 0 : msens-1 ] *phi; 

dsp= [ 1 : msens ] *dspace* (1 /msens ) ; 

%dsp=[l 2 5 11 13 16 18] *1/18; 
for nnn=l: msource 
• cira (nnn, : ) =ql ; 
dl (nnn, : ) =dsp; 

end 

for nb=l: msens 

cb ( : , nb ) =bear ing ’ ; 
end 

amat2= ( (2*pi*dl . *cos (cira-cb) ) ) ' ; 

.for ns=l: msource 

a 3 ( : , ns ) =amat2 ( : , ns) *sin (ang (ns ) *pi/18 0 ) ; 

end 

amat 3=exp ( j *a3 ) ; 

smat2=smat ^amatS . ' ; 

[nsamp, msens] = size(smat2); nsource=msource; 
namp=l/ (sqrt (2) *10 A (snr/20) ) ; 
for n=l : msens 

nmat ( : , n) =randn (nsamp, 1 ) + j *randn (nsamp, 1 ) ; 
end 



smat2=smat2+namp + nmat ; 

ymat=smat2;% estimate the spatial matrix 

cmat = zeros (msens, msens ) ; 

ymat = ymat - ones (nsamp, 1 ) + mean(ymat); % remove mean 

cmat = conj(ymat' + ymat) / nsamp; 

[umat, smat, vmat] = svd(cmat); 
svec = diag(smat); 
format compact 

% estimate bearing spectra 

prmat = zeros (msens , msens ) ; 
for i=l: msens 

if (svec(i) > 0) 

prmat = prmat + vmat ( : , i ) * vmat ( : , i ) ' / svec ( i ) ; 

end 

end 

prvec = prmat (:,1); 
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(pi/180) ; 



dtheta=l ; 

theta = [0 : dtheta: thr] ’ * 
va=[0: 1: 90] * (pi/180) ; 
pp= (length (theta) ) ; 
for nnl=l:pp 

ciral (nnl , : ) =ql ; 
d (nnl , : ) =dsp; 

end 

for nbl=l:msens 

cbl ( : , nbl ) -theta; 
end 

omega = 2*pi *d . *cos (ciral-cbl ) ; 
mth =’ length (theta) ; 

mlc = zeros (mth, 1 ) ; 

gvec = umat ( 1 , 1 : nsource ) . ' ; 
gmat = umat ( 2 : msens , 1 : nsource ) ; 

arvec = [1 ; - conj(gmat) * gvec / (1 - gvec'*gvec) ] 

tmat = vmat (:, nsource+1 : msens ) ; 
wvmat = tmat * diag (ones (msens-nsource, 1 ) . / 

svec (nsource+1 : msens ) ) * tmat*; 

tmat = tmat * tmat ’ ; 



nv=length ( va) ; 



% > search angle< % 

for bb=l:nv 
for k =l:mth 

steer = exp(j * omega (k, :)’ *sin (va (bb) )) ; 
beaml(k,bb) =1./ abs (steer* + inv( cmat) * steer); 
sminl(k,bb) = 1./ abs (arvec. ' * steer) /s 2; 
smusl(k,bb) = 1./ (steer* + tmat * steer); 
spisl(k,bb) = 1./ abs(steer* + vmat (:, msens) )^2; 
end 
end 



val= [-90 : 1 : 90 ] * (pi/180) ; 
mtl=length (val) ; 

% display the angular spectra 

theta = theta* 180/pi ; 
va=va* 180/pi ; 
val=val* 180/pi ; 

top=0 . 95 ; spacing=0 . 02 ; 
left-0-80; 
btnWid=0 . 1 ; 
btnHt=0 . 04 ; 
bottom=0 . 05; 
figure ( 1 ) 

uicontrol (* Units points ’ , ... 

' BackgroundColor ’ , [0 . 8 0.8 0.8], ... 

' FontSize *,19, ... 

' Fore groundColor *,[0.8 0.8 0.8], ... 

’ Position* , [178 264 100 19], ... 

* String * , * HOLD ON ! ! ' , ... 
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' Style ' , ' text 1 , ... 

' Tag ’ , ' StaticText2 ' ) ; 

fa=get(0, ' current figure ’ ) 
figure ( fa+1 ) , 

tname - [’SPIRAL CIRCULAR ANTENNA ARRAY,SNR=',int2str(snr)J; 

set ( gcf , ' name ' , tname ) 

asl=10*logl0 (abs ( (smusl + eps) ' ) ) ; 

as2 = 10*logl0 (abs ( (spisl+eps) ’ ) ) ; 

as3 = 10'*'logl0 (abs ( (sminl+eps) ’ ) ) ; 

as4=10* loglO (abs ( (beaml+eps) ’ ) ) ; 



bb=10*logl0 (abs ( (smusl + eps) ’ ) ) ; 
bias=min (min (bb) ) / 
as=max (max (bb) ) ; 

if thr==3 60 

mesh ( theta ’, va, bias+10*logl0 (abs ( (smusl+eps)’))); 

%title ([’ MUSIC, spiral circular array , SNR= ’ , num2str ( snr , 3 ) , * , max 

spacing : ’ , num2str (dspace, 3 ) ] ) ; 

axis([0 PO^thrv 0 90 -20 as]); 

zlabel ( ’power (dB) ’ ) ; 

xlabel ( ’ AZ (deg) ’ ) 

ylabel ( ’EL (deg) ’ ) 

else 

subplot (221) , mesh ( theta’ , va, bias + 10*logl0 (abs ( (smusl+eps) ’ ) ) ) 

% tit le ( ’ music ’ ) ; 

axis([0 90*thrv 0 90 -20 as]); 

bb=10*logl0 (abs ( (spisl + eps) * ) ) ; 

bias=min (min (bb) ) ; 

as=max (max (bb) ) ; 

zlabel ( ’ power (dB) ’ ) ; 

title ( ’MUSIC’ ) ; xlabel ( ’AZ (deg) ’ ) 

ylabel ( ’ EL (deg) ’ ) 

subplot ( 222 ), mesh ( theta ’, va, bias + 10*logl0 (abs ( (spisl + eps)’))) 
title ( ’ PHD’ ) ; 

zlabel ( ’ power (dB) ’ ) ; xlabel ( ’ AZ (deg) ’ ) ; ylabel { ’ EL (deg) ’ ) 
axis([0 90'*'thrv 0 90 -20 as] ) ; 

bb=10*logl0 (abs ( (sminl + eps) ’ ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

%xlabel ( ’Az-angle’ ) ; 

subplot ( 223 ), mesh ( theta ’, va , bias+ 10*logl0 (abs ( (sminl+eps)’))) 

title ( ’Mini-norm’ ) ; 

zlabel ( ’ power (dB) ’ ) ; 

title ( ’MUSIC’ ) ; xlabel ( ’ AZ (deg ) ’) 

ylabel ( ’EL(deg) ’ ) 

axis([0 90*thrv 0 90 -20 as]); 

%xlabel ([’ uniform circular array , SNR= int2str ( snr )]) ; 

% ylabel ( ’ El -angle ’ ) ; 
bb=10*logl0 (abs ( (beaml + eps) ’ ) ) ; 
bias=min (min (bb) ) ; 
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subplot (2X4 ; , mesh theta ’ , va , abs (bias) +10*logl0 (abs ( (beaml+eps, ' ) ; ) ; 
title ('Capon MLM ' ) ; 

2 label ( ' power (dB) ’ ) ; 
xlabel ( ' AZ (deg) ' ) 
ylabel ( 1 EL (deg) ' ) 
end 

AZ=be a ring* 180/pi ; 

EL=ang; 

title ( [ ' SPIRAL ANTENNA ‘ NUMBER= ' , int2str (msens) , ' SOURCES 
NUMBER= ' , int2str (msource ) ] ) 
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APPENDIX D 



SPECTRUM CALCULATION FOR C I RCULAR+L INEAR ARRAY 



dspacel=str2num (get (bdspacel, ’string’ } ) 
lm=str2num (get (bdl, ’ string ’ ) ) 
bearing=bearingn; 
msource=2 ; 
if nsource>6 
msource=6; 
else 

msource=nsource; 

end 

smat = zeros (nsamp,msource) / 

for n = lrmsource 

u = rand (ns amp, 1 ) -0.5; 
smat (:,n)= -sign(u) .* log ( 1-2 *abs (u ) ) ; 

end 

% the propagation matrix 

bearing = bearing * pi / 180; % degrees to radian 

msensm=msens 1 1 ; 

phi=2*pi/msensm; 

ql= [ 0 : msensm-1 ] *phi ; 

cira=ones (msource, 1 ) *ql ; cb=bearing ' *ones ( 1 , msensm) ; 

amat2= ( (2*pi*dspace*cos (cira-cb) ) ) ' ; 
amat22= ( (2*pi* [Or lm-1] ' + dspacel ) ) ; 

for ns=l:msource 

a 3 ( : , ns) =amat2 ( : , ns ) *sin (ang (ns ) * pi/180 ) ; 
a 32 ( : , ns ) =amat22*cos (ang (ns ) * pi /I 80 ) ; end 
a3= [a3 ; a32 ] ; 

amat3=exp ( j + a3 ) ; 

% the noisefree sensor signal 
smat2=smat *amat3 . ' ; 



% add white noise : 

[nsamp, msens] = size(smat2); 

nsource=msource; 

namp=l/ (sqrt (2) *10^ (snr/20) ) ; 
for n=l: msens 

nmat ( : , n) =randn (nsamp, 1 ) + j *radn (nsamp, 1 ) ; 
end 

nmat = nmat - ones (nsamp, 1 ) * mean(nmat); 

smat2=smat2 + namp'* r nmat ; 
ymat=smat2 ; 

% estimate the spatial cross-cumulant matrix 
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cmat = zeros (msens , msens ) ; 

ymat = yxnat - ones { nsamp, 1 , * mean (ymat;; remove mean 

cmat = conj(ymat' * ymat) / nsamp; % correlation matrix 

. determine number of sources 

[umat, smat, vmat] = svd(cmat); 
svec = diag(smat); 
s=[l 2 3] ; 
dtheta=l ; 
format compact 

estimate bearing spectra 

prmat = zeros (msens , msens ) ; 
for i=l:msens 

if (svec(i) > 0) 

prmat = prmat + vmat(:,i) * vmat(:,i)’ / svec(i); 

end 

end 

prvec = prmat (:, 1) ; 

theta = [0 : dtheta : thr ] ’ * (pi/180); 

va= [0:1:90] * (pi/180) ; 
pp= (length ( theta ) ) ; 

ciral=ones (pp, I) *ql; 

cbl=theta*ones (l,msensm) ; 

omega = 2*pi*dspace*cos (ciral-cbl ) ; %< 00000*^ 

mth = length (theta) ; 

mlc = zeros (mth, 1 ) ; 

gvec = umat ( 1 , 1 : nsource ) . f ; 
gmat = umat (2 rmsens, 1 rnsource) ; 

arvec = [1 ; - con j (gmat) * gvec / (1 - gvec'^gvec) ]; 
tmat = vmat (:, nsource + 1 rmsens ) ; 

wvmat = tmat * diag (ones (msens-nsource, 1 ) ./ svec (nsource+1 : msens ) ) * 

tmat ? ; 

tmat = tmat * tmat ' ; 
nv=length ( va ) ; 
if thr==3 60 
for bb=l:nv 
for k =l:mth 

steerl = exp(j * omega ( k, :)' *sin (va (bb ))) ; 
steer 2= exp ( j *cos (va (bb) ) *2*pi*dspacel* [0: lm-1] ' ) ; 
steer= [steerl ; steer2] ; smusl(k,bb) = 1./ (steer 1 * tmat * steer); 
end 
end 

else 

for bb=l : nv 
for k =l:mth 

steerl = exp(j * omega (k, :)’ *sin (va (bb) )) .; 
steer2= exp(j *cos (va (bb) ) *2 *pi *dspacel * [Orlm-l]'); 

steer= [steerl ; steer2] ; 
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beaml(k,bb = 1 . /abs ( steer ' + inv( cma t ) * steer) ; 
sminl(k,bb) = 1./ abs (arvec. ' * steer) A 2; 
smusl(k,bb) = 1./ (steer' * tmat + steer); 
spisl(k,bb) = 1./ abs (steer' + vmat ( : , msens ) ) A 2; 

end 
end 
end 

val=[-90: 1: 90] *(pi/180); 
mt l=length ( val ) ; 



o display the angular spectra 

theta = theta*180/pi; 

va=va* 180/pi ; 
vai=val * 180/pi; 



top=0 . 95 ; spacing=0 . 02 ; 
lef t=0 .80; 
btnWid=0 . 1 ; 
btnHt=0 . 04 ; 
bottom=0 . 05 ; 
figure ( 1 ) 

uicontrol ( ' Units ’ , ' points ' , ... 

' BackgroundColor ' , [0 . 8 0.8 0.8], ... 

' FontSize ',19, ... 

' ForegroundColor ' , [0 . 8 0.8 0.8], ... 

' Posit ion ' , [ 17 8 264 100 19], ... 

'String' , 'HOLD ON ! ! ' , ... 

' Style ' , ' text 1 , ... 

' Tag ' , ' StaticText2 ' ) ; 



fa=get(0, ' current figure ') ; 
figure ( f a+1 ) , 



asl=10*logl0 (abs ( (smusl+eps) 1 ) ) ; 
bb=10'* r loglO (abs ( (smusl + eps) * ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

if thr==360 

mesh( theta 1 , va, bias + 10*logl0 (abs ( (smusl + eps) '))); 

% title ( [ 1 MUSIC, linear+circular array, SNR= 1 , num2str ( snr, 3 ) , ' , C- 
s pacing : ’ , num2st r ( dspace , 3 ) , ’ , L- space : ' , num2str (dspacel , 3 ) , 1 , L-ant 
no : ' , nun\2str ( lm, 3 ) ] ) ; 
axis([0 90*thrv 0 90 -20 as]); 

%xlabel ( [ 1 ANTENNA NUMBER= int2str (msens ) ] ) ; 

%y label ( [ ' source number= ' , int2str (nsource ) ] ) ; 

%zlabel ( 'power (dB) ' ) ; 
z label ( ' power ( dB) ' ) ; 
xlabel ( ' AZ (deg) ’ ) 
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ylabel ' ' EL (deg) ' ) 
else 

as2=10*logl0 (abs ( ispisl+eps) ' ) ) ; 
as3=10 + logl0 (abs ( (sminl+eps) ' ) ) ; 
as4=10 + logl0 (abs ( (beaml+eps) ' ) ) ; 

subplot (221) , mesh ( theta 1 , va, bias+10 * loglO (abs , smusl+eps ) ' ; . ; 

axis ( [ 0 90*thrv 0 90 -20 as]); 

z label ( ' power (dB) ' ) ; 

title ( 'MUSIC' ) ; xlabel ( 'AZ (deg) ' ) 

ylabel ( ' EL (deg) ' ) 



bb=10*logl0 (abs ( (spisl+eps) ' ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

subplot (222 ) , mesh ( theta ',va,bias+10*logl0 (abs ( (spisl+eps )'))) ; 

zlabel ( ' power (dB) ' ) ; 

title ( ' PHD ' ) ; xlabel ( ' AZ (deg) ' ) 

ylabel ( 'EL(deg) ’ ) 

axis([0 90*thrv 0 90 -20 as] ) ; 
bb=10*logl0 (abs ( (sminl+eps) ' ) ) ; 
bias=min (min (bb) ) ; 
as=max (max (bb) ) ; 

subplot ( 223 ) , mesh ( theta ' , va, bias + 10 + log 10 (abs ( (sminl + eps )'))); 

title ( 'Mini-norm' ) ; 

zlabel ( ' power (dB) ' ) ; 

xlabel ( 'AZ (deg) ' ) 

ylabel ( ' EL (deg) ' ) 

axis([0 90*thrv 0 90 -20 as] ) ; 

tname= ( [ ' linear+circular , SNR= ' , int2str (snr ) ] ) ; 

%ylabel ( 'El (degree) ' ) ; 

bb=10*logl0 (abs ( (beaml+eps) ' ) ) ; 
bias=min (min (bb) ) ; 

%as=max (max (bb) ) ; 

subplot (22 4 ) , mesh ( theta ',va, abs (bias ) +10*logl0 (abs ( (beaml+eps ) ' ) ) ) 

title ('Capon MLM'); 

zlabel ( ' power (dB) ’ ) ; 

xlabel ( 'AZ (deg) ' ) 

ylabel ( 'EL (deg) ’ ) 

end 

AZ=bearing* 180 /pi ; 

EL=ang; 

tname= ( [ ' linear+circular , SNR=' , int2str (snr) ] ) ; 
set(gcf, ' name ' , tname ) 
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APPENDIX E 



SPECTRUM CALCULATION FOR LINEAR ARRAY 



clear nmat 

[nsamp, msens] = size(smatl); 
nsource=msource; 

namp=l/ (sqrt (2) *10^ (snr/20) ) ; 
for n=l : msens 

nmat (:,n) = randn (nsamp, 1 ) + j * randn ( nsamp, 1); 



end 

smat l=smat l+namp*nmat ; 
ymat=smat 1 ; 

1 estimate the spatial cross-cumulant matrix 

cmat = zeros (msens , msens) ; 

ymat = ymat - ones (nsamp, 1 ) * mean(ymat); % remove mean 

% determine number of sources 

cmat = conj(ymat' * ymat) / nsamp; % correlation matrix 



determine number of sources 

[umat, smat, vmat] = svd(cmat); 
svec = diag(smat); 



% estimate bearing spectra 

prmat = zeros (msens , msens ) ; 
for i=l:msens 

if (svec(i) > 0) 

prmat = prmat + vmat(:,i) * vmat(:,i)' / svec(i); 

end 

end 

prvec = prmat (:,1); dtheta=0.5; 



theta = [-90 : dtheta : 90 ] ' * (pi/180); 
omega = 2*pi *dspace*sin ( theta ) ; 



mth = length (theta) ; 

mlc = zeros (mth, 1 ) ; smus = mlc; spis = mlc; seig = mlc; par 
mlc ; 

beam = mlc; smin = mlc; 

gvec = umat ( 1 , 1 : nsource ) . ’ ; 
gmat = umat ( 2 : msens , 1 : nsource ) ; 

arvec = [1 ; - con j (gmat) * gvec / (1 - gvec ' *gvec ) ]; 

tmat = vmat ( : , nsource+1 :msens ) ; 
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wvmat tmat * diag (ones (msens-nsource, 1 j . / svec (nsource+1 : msens ) ) * 

tmat ' ; 

tmat = tmat * tmat 1 ; 
for k =l:mth 

steer = exp(j * omega (k) * [0 :msens-l] ' ) ; 

obeam(k) « abs {steer' * cmat * steer); 
beam{k) = 1 . /abs { steer ' * inv(cmat) * steer); 

smin(k) = 1./ abs(arvec. ' * steer) *2; 

smus(k) =1./ (steer* * tmat * steer) ; 
spis(k) = 1./ abs (steer ' * vmat {:, msens ) ) A 2; 

end 



o display the angular spectra 

theta = theta*180/pi; 

spec — abs ( [seig, smus, spis, mlc, par, smin, beam] ); 

% normalize to abs-max of unity for display only 

spmax = max(spec); 

spmax = ones (1,7) ./ spmax; 

spec = spec * diag(spmax); 

figure (2 ) 

Isubplot { 421 ) , stem(svec), title(tname) 

subplot ( 421 ) , plot (theta, 10*logl0 (spec ( : , 3) ) , ' r 1 ) 

title (' Uniform linear array 1 ) 

ylabel { ' power (dB) 1 ) 

gtext ( ' PHD' ) 

subplot ( 423 ) , plot (theta , 10* log 10 (spec ( : , 2 ) ) , ' r ' ) 

ylabel { ' power (dB) 1 ) , gtext ( 1 MUSIC ' ) 

subplot ( 4 25 ) , plot (theta , 10* log 10 ( spec ( : , 6) ) , ’ r ' ) 

ylabel ( ' power (dB) 1 ) , gtext ( ' Mini-norm 1 ) 

subplot { 427 ) , plot (theta, 10*logl0 (spec ( : , 7 ) ) , ’ r ’ ) , hold 
ylabel ( ' power (dB) ’ ) , gtext ( ' Capon MLM’ ) , xlabel ( ’ AOA (degree) ' ) 
tname = ['UNIFORM ANTENNA and NON-UNIFORM ANTENNA LINEAR ARRAY’] 
set(gcf, ' name 1 , tname ) 
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APPENDIX F 



- SPECTRUM CALCUIiATION FOR NON-UNIFORM LINEAR 
ARRAY 



clear nmat 

[nsamp, msens] = size(smat2); 
nsource=msource ; 

%snr=input { 1 snr= 1 ) ; namp=l/ ( sqrt ( 2 ) * 10 A { snr /20 ) ) ; for n=l: msens 
nmat ( : , n) =randn (nsamp, 1 ) + j *randn (nsamp, 1 ) ; 

end 

smat2=smat2+namp*nmat ; 
ymat=smat 2 ; 

% estimate the spatial cross-cumulant matrix 

cmat = zeros (msens , msens ) ; 

ymat = ymat - ones (nsamp, 1 ) * mean(ymat); % remove mean 

cmat = conj(ymat' * ymat) / nsamp; 

% correlation matrix 



% determine number of sources 

[umat, smat, vmat] = svd(cmat); 
svec - diag(smat); 
hold off 

% estimate bearing spectra 

prmat = zeros (msens , msens ) ; 
for i=l:msens 

if (svec(i) > 0) 

prmat = prmat + vmat(:,i) * vmat(:,i) f / svec(i); 

end 

end 

prvec = prmat (:,1); 



dtheta=0 . 5 ; 

theta = [ -90 : atheta : 90] T * (pi/180); 
omega = 2*pi*dspace*sin (theta ) ; 
mth = length (theta ) ; 

mlc = zeros (mth, 1 ) ; 

smus = mlc; spis = mlc;seig = mlc; par = mlc; 
beam = mlc; smin = mlc; 

gvec = umat (1, l:nsource) . ’ ; 
gmat = umat (2 :msens , 1 : nsource) ; 

ar.vec = [1 ; - conj(gmat) * gvec / (1 - gvec ' *gvec) ]; 

tmat = vmat ( : , nsource+1 rmsens ) ; 

wvmat = tmat * diag (ones (msens-nsource, 1 ) ./ svec (nsource+1 : msens ) ) * 

tmat f ; 

tmat = tmat * tmat f ; 
for k =l:mth 
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steer - exp(sqrt(-l) * omega(k) * asp'); 
,beam{k) = abs (steer' * cmat * steer); 
beam(k) = 1 . / abs (steer* *inv( cmat) ‘ steer); 
smin(k) - 1./ abs (arvec. ' * steer) A 2; 
smus(k) - 1./ (steer' * tmat 4 steer); 
spis(k) ~ 1 ./ abs (steer* * vmat ( : , msens ) ) A 2; 

end 



o display the angular spectra 

theta = theta*180/pi; 

spec = abs ( [smus, spis, smin, beam] ); 

: normalize to abs-max of unity for display only 

spmax = max (spec); 

spmax = ones (1,4) . / spmax; 

spec = spec * diag (spmax); 

figure ( 2 ) 

subplot { 4 22 ) , plot (theta, 10*logl0 ( spec ( : , 2) ) ) ; 

%set(gca, ’position*, [0.06 0.75 0.35 0.17]); 

title ( * Non-uniform linear array' ) , ylabel ( ’power (dB) ' ) , gtext ( * PHD* ) 
subplot (424 ) , 

plot (theta , 10* log 10 ( spec ( : , 1 ) ) ) , ylabel ( ’power (dB) * ) , gtext ( ' MUSIC ' ) 
subplot (426) , plot (theta, 10* log 10 ( spec ( : , 3 ) ) ) ; 
ylabel ( ' power (dB) ' ) , gtext ( *Mini-norm* ) 

%set (gca, ’ position *, [0 . 06 0.28 0.35 0.17]); 

% title ( 'Minimum norm method') 

subplot (428) , plot (theta, 10*logl0 (spec ( : , 4) ) ) ; 

ylabel ( 'power (dB) ' ) , xlabel ( ’AOA( degree) * ) , gtext ( 'Capon MLM ' ) 
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